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ABSTRACT 


Chaos  metrics  are  examined  as  a  tool  to  analyze 
atmospheric  three-dimensional  dispersion  models  at  the 
individual  particle  rather  than  the  aggregate  level.  These 
include  the  self-affine  fractal  dimension,  DA/  Shannon 
entropy,  S,  and  Lyapunov  exponent,  X.  Inter comparison  of  these 
metrics  is  first  performed  with  the  one-dimensional  logistics 
difference  and  the  two-dimensional  Henon  systems  of  equations. 
The  fractal  dimension  and  Shannon  entropy  are  then  measured  as 
a  function  of  the  inverse  Monin-Obukhov  length  (1/L)  for  two 
three-dimensional  Lagrangian  particle  dispersion  models,  the 
McNider  particle  dispersion  model  and  the  NPS  particle 
dispersion  model  now  under  development.  The  fractal  dimension 
and  Shannon  entropy  uncover  weaknesses  in  both  models  which 
are  not  obvious  with  standard  geophysical  measures.  They  also 
reveal  similarities  and  differences  between  the  atmospheric 
models  and  simple  chaos  systems.  Combined,  these  chaos 
measures  may  lend  detailed  insight  into  the  behavior  of 
Lagrangian  Monte  Carlo  dispersion  models  in  general. 
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I .  INTRODUCTION 


A.  OVERVIEW 

Unlike  the  relatively  simple  linear  relationships  studied 
in  basic  physics,  natural  phenomena  often  exhibit  complex 
nonlinear  behavior.  One  familiar  consequence  is  the  inability 
of  meteorologists  to  make  accurate  long-term  weather  forecasts 
(Devaney,  1990) .  Like  turbulence  in  other  systems,  atmospheric 
weather  systems  seem  to  follow  patterns  which  are  chaotic  and 
unpredictable  in  the  long  term.  The  behavior  of  a  particle 
released  in  such  a  system  is  deterministic  rather  than 
stochastic.  However,  small  differences  in  initial  conditions 
result  in  particle  paths  which  diverge  exponentially  with 
time.  Thus,  since  the  resolution  and  accuracy  of  measurements 
of  the  initial  conditions  are  always  limited,  forecasts  will 
diverge  from  reality  exponentially  with  time  such  that  long 
term  predictability  becomes  impossible  (Gleich,  1987) . 

In  fluids,  predictability  relates  to  particle  diffusion. 
Although  the  relationship  is  not  simple,  less  predictable 
systems  tend  to  more  diffusive. 

A  good  Lagrangian  particle  model  which  can  accurately 
simulate  diffusion  over  a  wide  range  of  atmospheric  conditions 
is  of  value  in  predicting  toxic  dispersion  from  various 
sources,  such  as  biochemical  warfare  agents,  Navy  LNG  tanker 


1 


leaks,  nuclear  reactor  and  weapons  plumes,  industrial  chemical 
plants,  large  rocket  and  missile  launch  emissions,  aborts,  and 
static  test  firings,  and  a  variety  of  planned  and  emergency 
releases  from  tanks  and  hoses  involved  in  storage,  transfer, 
and  transport  of  dispersive  volatile  toxins.  Presently, 
predictions  using  existing  particle  models  may  vary  greatly 
from  reality  when  operating  in  certain  atmospheric  conditions. 
(Venkatram  and  Wyngaard,  1988) . 

Heretofore,  atmospheric  particle  dispersion  models  have 
generally  been  compared  against  turbulence  data  in  terms  of 
spectra  and  standard  geophysical  turbulence  measures,  such  as 
the  turbulent  kinetic  energy  (tke) ,  its  component  velocity 
variances  (ctu2,  ct„2,  and  <xw2)  ,  and  the  Brunt-Vaisala  frequency 
( BVF) ,  an  atmospheric  stability  measure.  The  aggregate 
diffusion  of  many  Lagrangian  particles  has  also  been  compared 
with  the  results  of  Gaussian  plume  statistical  theories  as 
well  as  measurements  of  dispersing  clouds  (Venkatram  and 
Wyngaard,  1988) .  However,  complex  but  wholly  predictable 
periodic  behavior  such  as  ocean  tides  may  appear  quite  similar 
to  truly  chaotic  or  random  behavior  when  using  such 
techniques.  Another  potential  set  of  tools  in  analyzing 
particle  dispersion  model  performance  are  the  recently 
developed  chaos  metrics  such  as  fractal  dimension,  DA,  Shannon 
entropy,  S,  and  the  Lyapunov  exponent,  X  (Baker  and  Gollub, 
1990)  . 
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B.  OBJECTIVE 


This  study  examines  the  self-affine  fractal  dimension,  DA , 
and  compares  it  to  the  Shannon  entropy,  S,  and  the  Lyapunov 
exponent,  X.  These  metrics  are  first  applied  to  the  1-D 
chaotic  system  known  as  the  logistics  difference  relation  and 
the  chaotic  2-D  Henon  system.  These  simple  well  characterized 
systems  are  studied  to  determine  similarities  and  differences 
between  the  above  three  chaos  metrics.  Two  3-D  Monte  Carlo 
Lagrangian  scattering  routines  designed  to  simulate 
atmospheric  dispersion  are  then  studied  as  representative  of 
more  complex  real  world  systems. 

C.  WHY  USE  CHAOS  METRICS 

There  are  at  least  two  good  reasons  to  try  chaos  metrics 
on  dispersion.  One,  chaos  metrics  may  provide  a  means  to 
discriminate  between  periodic  wave  behavior,  chaos,  and 
differing  degrees  of  turbulence.  Second,  chaos  metrics  may 
provide  some  additional  insight  into  the  behavior  of  particle 
dispersion  in  3-D  scattering  routines. 

Kamada  and  DeCaria  (1992)  have  shown  that  though  nocturnal 
periodic  gravity  waves  have  quite  different  dispersion 
characteristics  from  turbulence,  the  two  cannot  be 
distinguished  readily  by  using  standard  atmospheric  turbulence 
measures  such  as  the  Brunt-Vaisala  frequency  (BVF) ,  vertical 
velocity  variance  (aw2)  ,  turbulent  kinetic  energy  (tke)  , 
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buoyancy  length  scale  (lb),  or  spectral  analysis  using  FFTs. 
The  self-affine  fractal  dimension,  DA ,  was  shown  to  be  the 
only  facile  wave/turbulence  discriminant  tested.  Since  the 
Shannon  entropy  and  Lyapunov  exponent  are  two  other  standard 
chaos  measures,  they  might  also  be  tested  as  possibly  useful 
dispersion  measures. 

In  the  3-D  Monte  Carlo  scattering  routines  which  model 
atmospheric  dispersion,  such  as  the  McNider  Particle 
Dispersion  Model,  chaos  metrics  might  quickly  determine 
certain  model  characteristics  for  a  given  range  of  parameters. 
For  example,  in  examining  an  expanding  cloud  of  particles,  the 
Lyapunov  exponent  would  measure  the  divergence  rate  of  the 
particles.  The  Shannon  entropy  would  measure  the  evenness  of 
distribution  of  particles  within  the  expanding  cloud.  The 
self-affine  fractal  dimension  could  measure  the  total  apparent 
distance  a  particle  travels  in  a  time,  T,  as  a  function  of  the 
time  resolution,  e,  used  within  T.  This  might  give  an 
indication  of  the  range  of  scales  over  which  mixing  occurs. 
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II.  THEORY 


A.  OVERVIEW 

There  are  several  methods  of  measuring  chaos  and 
turbulence.  Of  interest  for  the  present  purposes  are  the  self- 
affine  fractal  dimension,  DA/  the  Shannon  entropy,  S,  and  the 
Lyapunov  exponent,  X.  Some  standard  geophysical  measures  will 
also  be  used  for  comparison  purposes  in  the  3-D  case. 
Descriptions  of  each  are  given  below  and  the  logistics 
difference,  Henon,  and  atmospheric  particle  diffusion  systems 
to  which  they  are  applied  are  then  described. 

B.  FRACTAL  DIMENSION  <DA) 

The  self-similar  fractal  dimension,  DB,  can  be  described 
readily  with  the  Cantor  set  (Devaney,  1990)  .  To  form  this  set, 
a  straight  line  is  drawn  with  the  middle  third  removed,  as  in 
the  second  line  of  Figure  1.  The  two  resulting  straight  lines, 
which  are  one-third  the  original  line  length,  are  then 
similarly  subdivided.  The  four  new  lines,  all  l/9th  the 
original  line  length,  are  further  subdivided,  ad  infinitum . 
From  top  to  bottom  in  Figure  1,  the  resulting  sets  of  lines 
are  self-similar ;  that  is,  the  manner  in  which  the  geometry  is 
altered  is  repeated  at  all  successive  levels  of  resolution. 

The  Koch  snowflake  is  another  geometrically  self-similar 
figure  (Devaney,  1990).  It  is  constructed  by  removing  the 
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Figure  1 ,  Construction  of  the  Cantor  Set 
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middle  third  of  each  side  of  an  equilateral  triangle,  and 
replacing  it  with  two  pieces  of  equal  length,  creating  a  six- 
pointed  star.  This  star  has  12  sides,  all  of  length  1/3.  The 
middle  third  of  each  of  the  12  sides  is  again  removed,  and 
replaced  with  two  lines  of  length  1/9.  Again,  the  process  is 
repeated  ad  infinitum. 

Like  the  Cantor  set,  the  Koch  snowflake  is  also  self¬ 
similar  [Figure  2].  The  jagged  sides  appear  geometrically 
similar  at  increasing  levels  of  resolution. 

The  self-similar  fractal  dimension  is  defined  by  (Devaney, 
1990)  as 

D  _  In  ( Total  length  of  pieces)  (II-l) 

B  In  ( resolution ) 

In  the  Cantor  set  example,  the  total  length  of  segments  of 
unit  length  at  level  n  is  2",  and  the  resolution  level  is  3n, 
so  that 


D  =  ln  2”  =  n  |n---  -  0.6309  .  ( II-2 ) 

8  In  3"  n  In  3 

Similarly  for  the  Koch  snowflake,  there  are  4  pieces  for 
each  level  of  n  with  a  magnification  of  3,  so  that 

D  =  In  4n  a  12S2  .  •  ( 1 1  —  3 ) 

B  In  3n 

With  geometric  figures,  D0  is  unitless;  the  measure 
involves  a  length  divided  by  a  length.  However  this  definition 
of  fractal  dimension  cannot  apply  directly  to  a  time  series 
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(1990)  . 
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trace,  since  it  would  involve  the  square  root  of  the  amplitude 
squared  plus  the  time  squared,  i.e.. 


(II-4) 


This  definition  must  be  adjusted  when  applied  to  time  series 
data  to  ensure  that  the  end  result  is  independent  of  the 
arbitrary  unit  scaling  between  amplitude  and  time. 

There  are  other  ways  to  characterize  the  fractal  dimension 
of  a  system.  For  a  more  suitable  times  series  measure,  McHardy 
and  Czerny  (1987)  redefined  the  length  metric  as 


T 

Lie)  =  -f  lF(t+e)  -  Fit)  |  dt  .  (H-5) 

c  j 
0 

where  F  is  the  amplitude  of  the  time  series  at  time,  t,  e  is 
the  time  increment,  and  T  is  the  time  window  over  which  L  is 
defined . 

This  definition  effectively  avoids  the  units  scaling 
problem.  Since  the  inverse  time  in  1/e  cancels  the  time  units 
in  the  integral,  L(e)  is  only  in  units  of  amplitude.  The 
fractal  dimension  is  then  defined  as 


DA  =  -  d  j~n  L(e)  .  (II-6) 
d  In  e 

McHardy  and  Czerny  applied  these  definitions  to  the  time 
variance  of  X-ray  luminosity  data  from  the  Seyfert  galaxy 
NG5506.  Collins  and  Kastogi  (1989)  later  applied  McHardy  and 
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Czerny's  definitions  in  analyzing  gravity  wave  spectra  from 
the  atmospheric  mid-troposphere.  Recently,  Kamada  and  DeCar ia 
(1992)  applied  DA  as  a  tool  for  discriminating  waves  from 
turbulence  in  nocturnal  atmospheric  boundary  layer  data. 

To  actually  compute  this  function,  the  integral  is 

approximated  with  a  numerical  summation,  so  that 

Lie)  =  £|F(t-e)  -  F(t)|  .  (H-7) 

e  o 

Since  at  each  resolution  <St=e,  the  leading  term  on  the  right- 
hand  side  is  always  unity,  so  that  the  length  is  now 

T 

Lie)  =  £|F(t-e)  -  F(  t) |  .  (II-8) 

0 

Thus  L ( e )  is  the  total  amplitude  change  over  a  time  series  of 
length,  T,  for  a  given  time  resolution,  e.  In  this  form  it  is 
quite  clear  that  time  is  removed  from  the  length 

determination,  so  that  L(e)  does  not  depend  on  some  arbitrary 
scaling  between  amplitude  and  time. 

C.  SHANNON  ENTROPY  (S) 

1.  Definition 

When  a  particle  is  first  released,  its  initial 

location  is  known  and  completely  specified,  so  the  information 
entropy  (defined  later)  for  its  location  is  zero.  Later, 

according  to  given  equations  of  motion,  its  position  diverges 
from  the  initial  point.  If  the  range  of  its  possible  positions 
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is  partitioned  into  equal  segments,  then  for  a  given  time 
interval,  T,  its  motion  can  be  recorded  in  terms  of  occupancy 
time  for  each  segment.  Then  the  particle  dispersion  rate  might 
be  measured  by  the  seeming  degree  of  randomness  of  occupancy 
time  or  evenness  of  the  state  probability  distribution  over 
time  period,  T.  Shannon  or  information  entropy  is  defined  by 
the  state  probability  distribution,  so  entropy  may  be  regarded 
as  a  measure  of  dispersion  rate  for  a  time  series  of  particle 
states.  This  can  apply  to  the  actual  particle  position,  its 
velocity,  or  its  phase  velocity.  Shannon  entropy  is  defined 
simply  as  (Baker  and  Gollub,  1990) 

S  =  -  £  Pi  ln  Pi  #  (H-9) 

i»l 


where 


5  s  system  entropy, 

N  s  number  of  permitted  states,  and 

N 

pi  =  probability  of  state  i,  such  that  ^  pi  =  1.0 

i  =1 


Then,  for  N  permitted  states  (or  position  intervals) , 
the  maximum  possible  entropy  corresponds  to  equal  occupancy 
time  in  each  of  the  N  states,  so  Pj  is  1/N  for  all  i.  That  is. 


N 


■-E 


i  *1 


( II-10a) 


Expanding  the  summation  results  in 
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( II-10b) 


•^max  =  —  In—  +  —  In—  +  —  In  — 
N  N  N  N  N  N 


and  collapsing  the  common  multiples  gives 


5max  ‘  ~N  (  N  ln  n) 


^max  ~  ^  / 


(II-lOc) 


which  corresponds  to  total  randomness,  or  a  completely  even 
particle  distribution  across  all  allowable  states.  Also, 


*^min  ^  > 


(II-H) 


which  corresponds  to  all  particles  being  in  one  state. 

For  a  randomly  moving  particle  in  2-D  space,  the 
computation  is  similar.  The  2-D  space  may  be  divided  into, 
say,  a  100  X  100  grid.  If  the  particle  is  moving  completely 
randomly,  at  time,  t,  it  may  be  in  any  one  of  the  grid  boxes 
with  equal  probability.  The  entropy  for  this  system  is  then 


S  =  - 


— —  ln  — — )  +  ( — —  ln  — — 

)o2  1002/  v ioo2  loo2) 

S  =  -  1002  ln  — —  1  , 

i  ioo2  ioo2) 


which  is  the  same  as 


( II-12a) 


s  =  In  1002  =  SL 


(II-12b) 
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2 .  Shannon  Entropy  for  an  N-D  System 

For  a  1-D  system,  the  domain  is  simply  partitioned 
into  n  intervals,  and  the  probability  computed  for  each 
interval.  For  a  simple  system  such  as  recursion  of  the 
logistics  difference  equation, 

=  \ixn(l-xn)  ,  (11-13) 

the  probability  for  the  ith  interval,  p;,  is  simply  the  number 
of  times  the  interval  has  been  occupied  after  n  number  of 
recursion  steps,  divided  by  the  total  number  of  steps. 

Note  that  the  number  of  partitions  should  be 
appropriate  for  the  total  number  of  steps.  Having  too  few 
intervals  is  equivalent  to  too  low  a  resolution  of  the  entropy 
and  may  result  in  a  misleadingly  low  value  of  S.  For  instance, 
with  only  one  interval,  p;  =  1,  and  S  =  0.  Ideally,  the  number 
of  intervals  for  maximum  resolution  of  the  entropy  of  the 
system  is  ew,  where  X  is  the  positive  Lyapunov  exponent  for 
the  system  and  N  is  the  number  of  steps  (Baker  and  Gollub, 
1990,  pp. 126-129)  .  Since  X  is  typically  of  order  unity,  ew  can 
easily  be  a  computationally  intractable  number. 

Again,  for  3-D  systems,  assuming  N  equal  segments 
along  each  axis,  the  number  of  partitions  will  be  N3,  so 
computation  quickly  becomes  unwieldy  for  large  N.  In  practice, 
N  of  order  102  to  103  has  been  used  by  some  authors  (Baker  and 
Gollub,  1990,  pg.  88)  as  a  compromise  between  entropy 
resolution  and  computational  efficiency.  In  analyzing  real 
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data,  the  number  of  partitions  should  be  related  in  some 
direct  fashion  to  the  maximum  resolution,  e;,  of  the  measuring 
devices,  and  the  total  measurement  time,  T.  The  total  number 
of  partitions  should  probably  not  exceed  T/€u  or  otherwise 
even  a  completely  random  distribution  would  still  result  in 
some  intervals  unoccupied. 

3.  Ln  vs  I»og2 

Wolf  maintains  that  the  Shannon  entropy  should  be 
computed  with  log2  rather  than  the  natural  logarithm,  logc 
(Wolf,  1986,  pg.  276).  With  log2,  the  Shannon  entropy  relates 
directly  to  information  in  the  form  of  binary  bits,  since  one 
bit  of  information  can  be  specified  as  being  in  only  one  of 
two  states,  e.g.,  true  or  false,  on  or  off,  one  or  zero.  That 
is,  the  log2  basis  sets  the  entropy  equal  to  the  minimal 
length  of  binary  code  required  to  describe  the  state  of  the 
system.  If  all  particles  occupy  only  one  state,  turning  the 
bit  for  that  state  to  the  "on”  position  is  sufficient  to 
specify  the  state  of  the  system.  If  the  particles  are  evenly 
distributed  among  all  states,  the  length  of  binary  code 
required  to  specify  the  system  is  equal  to  the  number  of 
states,  which  means  that  the  system  is  completely  random 
(Tribus  and  Mclrvine,  1971).  However,  for  the  purposes  of  this 
study,  the  Shannon  entropy  can  also  be  written  in  the 
following  form: 
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(11-14) 


S  =  -kJ2  PilnPi 

i-l 

The  only  difference  in  measuring  entropy  by  the 
natural  logarithm  versus  base  two  is  the  value  imposed  on  the 
constant  K.  Most  computations  assume  that  K=l.  As  long  as  the 
method  of  computing  S  is  similar  when  making  comparisons  of 
computed  entropy,  there  need  be  no  confusion. 

D.  LYAPUNOV  EXPONENT  (X) 

For  an  expanding  turbulent  cloud  of  particles,  one  measure 
to  describe  the  system  is  the  particle  divergence  rate.  Due  to 
turbulence,  the  trajectories  of  two  particles  arbitrarily 
close  will  diverge  with  time.  The  divergence  rate  can  be 
characterized  by  Lyapunov  exponents.  The  Lyapunov  exponent,  X, 
is  also  a  measure  of  the  system's  sensitivity  to  initial 
conditions.  In  an  N-dimensional  system,  there  will  be  N 
Lyapunov  exponents;  however,  these  do  not  necessarily 
correspond  to  coordinate  axes.  So  in  a  Cartesian  coordinate 
system,  X,,  X2,  and  X3  are  locally  defined  Lyapunov  exponents 
which  generally  cannot  be  described  as  X„,  Xy,  and  Xz. 
Direction  is  adjusted  for  each  point  along  the  trajectory.  In 
one-dimension,  the  Lyapunov  exponent  is  defined  by  (Wolf, 
1986,  p.  275)  as 
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(11-15) 


X  =  lim  1  T  In  I f'U)  !  , 

AT--  AT  ffQ 

which  can  be  described  over  a  map  domain  as  an  integral, 

N 

X  =  f  Pi  In  |f'(i)  I  di  .  (11-16) 

0 

Overall,  the  Lyapunov  exponent  can  be  viewed  as  a  measure  of 
the  average  local  stretching  rate  of  the  particle  trajectory, 
as  indicated  by  the  log  of  the  length  of  the  slope,  |f'(i)|, 
weighted  by  the  probability  of  encountering  that  slope. 
(Wolf,  1986) 

The  Lyapunov  exponent  is  related  to  the  entropy,  S,  as 
well  as  the  information  loss  rate.  For  instance,  if  a  measured 
position  is  presently  known  to  a  precision  of  16  bits,  and 
A  =  2  bits  per  second,  then  the  particle's  future  trajectory 
cannot  be  predicted  to  any  degree  of  precision  beyond  8 
seconds  (16  orbits). 

It  should  be  noted  that  the  Lyapunov  exponent  defines  the 
average  rate  of  loss  of  predictive  power,  and  may  vary  locally 
along  the  orbit. 

The  Lyapunov  exponent  is  readily  computed  for  one 
dimension,  but  in  higher  dimensions  calculation  becomes  more 
complex.  In  three  dimensions,  the  directions  of  trajectory 
divergence  and  contraction  must  be  defined  in  terms  of  local 
tangents,  which  vary  from  point  to  point,  so  the  calculation 
of  the  exponents  must  constantly  adjust  for  each  change  in 
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direction.  Wolf  has  developed  algorithms  for  computing  higher 
dimensional  Lyapunov  exponents  which  involve  reorienting  the 
major  axis  of  an  ellipse  for  each  point,  then  renormalizing 
the  function  after  every  few  points  so  that  the  unit  ellipsoid 
does  not  overlap  an  attractor  (Wolf,  1986) . 

In  3-D,  three  Lyapunov  exponents  are  required  to  classify 
the  system.  A  negative  exponent  indicates  a  dissipative 
dimension.  If  all  three  exponents  are  negative,  the  system  is 
dissipative,  e.g.,  a  pendulum  settling  down  to  a  fixed  point. 
If  an  exponent  is  zero,  the  system  orbits  about  a  fixed  point. 
If  one  of  the  three  exponents  is  positive,  the  system  is 
chaotic;  the  orbital  trajectories  are  diverging. 

E.  GEOPHYSICAL  MEASURES  OF  TURBULENCE 

1.  Overview 

Geophysicists  use  several  standard  methods  to  describe 
turbulence.  Discussed  here  are  those  methods  used  to  examine 
the  McNider  and  NPS  dispersion  models. 

2.  Turbulent  Kinetic  Energy  (tke) 

The  mean  turbulent  kinetic  energy,  or  TKE,  is  defined 
as 


tke 


(11-17) 


where 

alt  s  standard  deviation,  u  component  of  particle  velocity, 
av  s  standard  deviation,  v  component  of  particle  velocity, 
a.M  =  standard  deviation,  w  component  of  particle  velocity. 
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3.  Gradient  Richardson  Number  (Ri) 

The  gradient  Richardson  number  is  defined  by 


(11-18) 


where  the  overbars  signify  time  averages,  and 

g  =  earth's  gravitational  acceleration, 
z  =  vertical  position  above  the  surface, 

6  =  potential  temperature,  i.e.  the  temperature 

normalized  for  adiabatic  expansion  with  pressure,  so 

/  \  ^ 

e  =  r  c*  ,  (ii-i9) 

\Po) 

where  Rd  is  the  gas  constant  for  dry  air,  and  Cp  is  the 
specific  heat  at  constant  pressure. 

Ri  is  used  in  place  of  the  Reynolds  number  as  a 
dynamic  indicator  of  turbulence  when  the  atmosphere  displays 
a  non-neutral  density  profile.  The  numerator  is  proportional 
to  the  buoyant  production  or  destruction  rate  of  tke, 
depending  on  negative  or  positive  sign,  respectively.  The 
denominator  is  proportional  to  the  shear  generation  rate  of 
tke,  which  is  nearly  always  positive.  Thus,  a  positive 
Richardson  number  indicates  a  stable  atmosphere  (increasing 
potential  temperature  with  height)  and  suppression  of 
turbulence.  Commonly,  when  the  buoyancy  destruction  rate  of 
tke  exceeds  1/4  the  shear  production  rate,  turbulence  is 
suppressed,  i.e.,  when  Ri  >  1/4  (Stull,  1988,  pg.  176). 
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4.  Brunt-Vaisala  Frequency  (BVF) 


The  Brunt-Vaisala  frequency,  BVF,  is  a  measure  of 
static  stability.  BVF  is  defined  by 


BVF  =  ,  (11-20) 

N  0v Sz 

and  in  principal  gives  the  highest  gravity  wave  frequency 
which  a  fluid  can  support.  It  is  undefined  for  negative 
temperature  gradients,  i.e.,  an  unstable  atmosphere. 
(Sorbjan,  1989,  pg.  35). 

5.  Buoyancy  Length  <1B) 

The  buoyancy  length,  lb,  is  the  standard  deviation  of 
the  vertical  velocity  divided  by  the  Brunt-Vaisala  frequency, 


7  s  w 
B  BVF 


(11-21) 


The  buoyancy  length  is  meant  to  be  a  measure  of  the  dominant 
eddy  scale.  (Stull,  1988,  pg.  310). 

6.  Monin-Obukhov  Length  (L) 

The  Monin-Obukhov  length  is  given  by 


ul  0 

L  =  - - - -  ,  (11-22) 

gkw'e  o 


where 
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o  s  surface  temperature  flux  , 
k  =  von  Karman  constant  =  0.4  , 


and  u. 


V 


(11-23) 


u.  is  the  friction  velocity,  a  measure  of  surface  drag 
due  to  turbulent  friction.  The  Businger  function,  \pm,  is  a 
stability  correction  which  is  approximated  with  a  rational 
function  by  (Kamada,  1992b) 


—  >  0 ( stable ) , 
L 


Vm  =  { 


0.15961583  -  5.4151107  — 

_ _  L 

1  -  3.59002232  -^  -  0.799168457j 


i) 


—  £  0  {unstable) 

jL / 


(11-24) 

In  the  convective  boundary  layer,  L  is  proportional  to 
the  height  at  which  the  buoyant  generation  rate  of  tke  matches 
the  shear  generation  rate.  L  is  the  primary  stability  measure 
for  the  atmospheric  surface  layer  and  is  also  used  in 
combination  with  the  inversion  height,  zx,  to  characterize 
boundary  layer  stability.  L  >  0  indicates  a  stable  boundary 
layer  where  vertical  turbulence  is  suppressed  by  positive 
density  gradients  with  height,  z.  L  <  0  indicates  an  unstable 
surface  layer  where  upward  vertical  motion  is  encouraged  by 
negative  density  gradients.  A  stable  atmosphere  implies  that 
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a  vertically  displaced  air  parcel  will  tend  to  return  to  its 
original  height.  (Businger,  1973) . 

F.  1— D  AND  2-D  CHAOS  EQUATIONS.  LOGISTICS  &  HENON  FUNCTIONS 
1.  Overview 

The  initial  focus  of  this  study  is  to  compare  and  test 
some  simple  potential  dispersion  metrics.  The  fractal 
dimension,  DA ,  the  Shannon  entropy,  S,  and  the  Lyapunov 
exponent,  X,  are  simple  expressions  which  can  monitor 
transitions  between  periodic  and  chaotic  behavior,  akin  to  the 
transition  between  laminar  and  turbulent  behavior  in  a  real 
fluid.  Two  such  expressions,  well  characterized  in  the 
literature,  are  the  logistics  difference  equation: 

xnn  =  iixn(l-xn)  ,  (11-25) 


and  the  Henon  system: 


xn+1  =  1  -  ax2n  +  yn  , 

Vn*  1  =  bxn  • 


(11-26) 


(Gould  and  Tobochnik,  1988,  pp. 152-178;  Baker  and  Gollub, 
1990)  . 

2.  The  Logistics  Difference  Equation 

Though  simple  and  one  dimensional,  recursive  iteration 
of  the  logistics  difference  equation  results  in  chaotic 
behavior  for  certain  parameters.  For  0  <  x0  <  1  and  n  <  1,  xn 
converges  to  0.  For  1  <  m  <  3 ,  and  the  same  range  for  x„,  xn 
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converges  to  m/4.  For  3  <  m  <  3.6,  x,,  fluctuates  between  2" 
discrete  points;  while  beyond  m  =  3.6,  the  solution  is  nearly 
random,  i.e.,  chaotic.  [Figure  3] 

3.  Graphical  Analysis  of  the  Logistics  Difference 

Equation 

The  logistics  difference  equation  is  parabolic,  so 
that  for  an  initial  x0  (given  as  0.04  with  m=2 . 9  in  Figure  4), 
drawing  a  vertical  line  to  the  parabola  corresponds  to  x,. 
Then  a  horizontal  line  drawn  to  the  inscribed  straight  line  of 
slope  unity  corresponds  to  recursing  x^,  to  a  new  value  for 
x„.  Reiterating  this  procedure  marches  the  solutions  to  a 
final  value  of  0.655. 

This  final  value  can  be  determined  exactly  from  the 
original  function.  At  steady  state,  xn+,  =  x„.  Then  for  Figure 
4 ,  where  m  =  2.9, 

=  2,9xn(l-xa)  ,  (11-25) 

and  solving  for  xn: 

xn  =  0 ,  X,  =  0.65517  .  (11-26) 

Note  since  the  equation  is  quadratic,  that  there  are 
two  solutions.  The  xn=0  solution  is  not  stable  (as  defined 
below).  The  general  rule  is:  if  the  magnitude  of  the  slope  of 
the  function  in  the  region  of  the  solution  is  greater  than 
45°,  the  fixed  point  is  an  unstable  solution;  if  the  slope  is 
less  than  45°,  the  fixed  point  is  a  stable  solution. 
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When  the  function  is  chaotic,  there  are  still  only  two 
solutions  to  the  quadratic  equation.  However,  both  points  are 
unstable:  they  both  repel  rather  than  attract  the  value  of 
xn+1.  The  two  points  can  still  be  computed,  resulting  in 

xn  =  0,  xn  =  0.74  .  (H-27) 

Another  way  of  looking  at  this  is  by  taking  the  derivative 
of  x,,.,,  with  respect  to  x„: 

df—  =  H  (1  -  2 Xn)  .  (11-28) 

d  xn 

Checking  the  chaotic  system  of  Figure  5,  when  xn  =  0.74,  the 
slope  is 

Wv 

:-n^  =  3.9  (1  -  2  (0.74))  -  -1.87  ,  (11-29) 

which  is  less  than  -1,  and  hence  unstable. 

This  concept  is  vital  when  visualizing  what  the  Lyapunov 
exponent  is  measuring.  Any  time  both  solutions  are  at  points 
where  the  slope  is  greater  than  45°  or  less  than  -45°,  the 
Lyapunov  exponent  is  greater  than  zero  (In  of  the  slope,  which 
is  greater  than  unity) ,  and  the  system  is  chaotic. 

4.  Bifurcation  Diagrams 

A  map  of  the  possible  values  of  x„  for  various  m  is 
called  a  bifurcation  diagram  (Figure  6].  With  the  bifurcation 
diagram,  the  complete  behavior  of  the  function  can  be 
determined.  For  example,  for  0  <  n  <  3,  xn  tends  to  only  one 


26 


/X 


Figure  6.  Logistics  Difference  Equation  Bifurcation  Diagram. 
Xn+I  =  MXn(l-xn)  . 
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value;  the  repetitive  recursion  of  xn+1  to  xn  points  to  one 
attracting  fixed  point,  corresponding  to  a  dissipative  system. 
At  m  >  3 ,  the  bifurcation  diagram  splits  into  two  fixed 
points.  The  solution  of  xn+1  bounces  back  and  forth  between 
these  two  points.  This  is  period  doubling  to  period  2.  Note 
for  still  larger  values  of  m  that  there  is  another  split  to 
period  4,  then  period  8,  followed  by  a  rapid  series  of 
bifurcations  culminating  in  chaos,  where  the  points  appear  to 
be  distributed  over  a  wide  range  of  xn. 

Another  curiosity  appears  in  the  bifurcation  diagram. 
There  are  numerous  "windows"  within  the  chaotic  region.  The 
most  obvious  is  around  m  =  3.82,  where  the  period  4  and  8 
behavior  seems  to  appear  again.  Increased  resolution  would 
show  many  more  such  windows.  Also  note  that  period  widths  get 
progressively  narrower,  and  that  for  higher  periods  it  is 
difficult  to  distinguish  say  period  64  from  chaos. 

5.  The  Henon  Equations 

As  with  the  logistics  difference  equation,  the  Henon 
set  has  stable  solutions  for  a  range  of  parameters  a  and  b 
which  correspond  to  one  fixed  point.  Other  values  for  a  and  b 
result  in  two  fixed  points,  or  period  doubling  [Figure  7],  and 
period  4  and  higher.  For  still  higher  values  of  a  and  b,  the 
result  is  chaotic  behavior:  the  possible  positions  appear  to 
be  spread  throughout  a  discrete  range  of  x  and  y  [Figure  8]. 


28 


The  Henon  bifurcation  diagram  resembles  the  logistics 
difference  equation  diagram,  but  has  some  obvious  differences 
[Figure  9].  One,  it  does  not  show  symmetric  splitting,  but 
rather  is  skewed  downward  from  the  centerline.  Two,  at  the 
value,  a  =  1.08  (with  b  =  0.3),  a  window  appears  in  the 
bifurcation  map  with  new  values  for  xn;  the  periodicity  in 
this  region  occurs  at  points  outside  the  previous  range  of  xn. 
These  new  points  seem  to  appear  out  of  nowhere,  and  are  not 
linked  to  previous  points  by  bifurcation. 

G.  3-D  ATMOSPHERIC  PARTICLE  DISPERSION  MODELS 
1.  Overview 

The  second  focus  of  this  study  is  to  examine  chaos 
metrics  as  performance  measures  for  3-D  Monte  Carlo  scattering 
routines.  Of  immediate  interest  is  particle  dispersion  within 
the  atmosphere.  The  McNider  Dispersion  Model  is  commonly  used 
to  estimate  the  transport  and  diffusion  of  atmospheric 
pollutants.  The  model  estimates  diffusion  based  on  atmospheric 
flow  predictions,  and  simulates  the  behavior  of  pollutant 
clouds  by  representing  them  as  the  ensemble  trajectories  of 
numerous  Lagrangian  point  particles  (Pielke,  1984) .  For  this 
study,  the  entire  range  of  possible  atmospheric  stabilities 
was  studied,  from  -0.2  <  1/L  0.1.  Therefore,  to  render  the 
computations  tractable,  the  flow  predictions  were  obtained 
from  boundary  similarity  theory  developed  by  Sorbjan  (1990) 
and  Kamada  (1992  a,b)  rather  than  from  a  mesoscale  windflow 
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model,  since  steady  state  mean  flows  could  be  assumed.  The 
difference  is  not  large,  since  most  of  the  similarity 
algorithms  are  also  used  in  the  turbulence  closure  scheme  for 
the  windflow  model. 

Several  weaknesses  were  perceived  in  the  McNider 
formulation.  Therefore,  a  similar  model  was  developed  at  NPS, 
and  was  also  tested.  It  features  a  double  Gaussian  vertical 
velocity  skewness  scheme  and  damping  of  the  particle 
oscillations  within  the  boundary  layer  (Kamada  1992b) . 

2.  The  McNider  Dispersion  Model 

For  this  study,  mean  flow  components  are  neglected  in 
order  to  focus  on  the  turbulent  fluctuations.  The  McNider 
Dispersion  Model  utilizes  the  form 

x(t  +  At)  =  x(t)  +  u(t)At, 

y  ( t  +  At)  =y(t)  +  v(  t)  A  t,  (11-30) 

z(t  +  At)  =  z(C)  +  w(t)At, 

where  x,  y,  and  z  define  the  particle  position  in  a  Cartesian 
coordinate  system,  and  the  u,  v,  and  w  terms  refer  to  the 
fluctuating  turbulent  velocity  components,  respectively.  The 
velocity  components  are  computed  by  (Pielke,  1984) 

u(t)  =  u{t  -  At)i?u(At)  +  u'(t  -  At), 

v(t)  =  v(t  -  At)£v(At)  +  v'(t  -  At)  ,  (11-31) 

w(  t)  =  w(C  -  At)i?„(At)  +  v/(t  -  At), 
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where  R  refers  to  the  Lagrangian  autocorrelation  factors  for 
each  velocity  component  and  is  a  function  of  At.  The  second 
terms  on  the  right  refer  to  random  fluctuating  components  of 
the  velocity,  and  are  related  to  the  flow  field's  turbulent 
kinetic  energy.  This  random  fluctuation  of  the  velocity  is 
constrained  by  the  fluid's  turbulent  kinetic  energy  with 

o'u  =  ouJl  -  i?u(At)  , 

a'v  =  ov\Jl  -  i?v( A t)  ,  (11-32) 

o'w  =  a w\j  1  -  Rl(bt)  , 


where  the  a  terms  refer  to  the  standard  deviations  of  the 
velocity  components.  These  parameters  are  not  directly 
computed  from  atmospheric  windflow  in  McNider's  model.  Instead 
an  empirical  formula  is  used  to  deduce  these  values: 


where 


(11-33) 


fO  .  35z,  z  <  205/n 

[^7 0 777,  Z  2  205m 


(11-34) 
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K,,,  is  the  local  exchange  coefficient,  described  later  in  the 
NPS  model  discussion.  To  account  for  advection  across  vertical 
gradients  of  crw. 


do. 


=  il:-5iwU  -  At>  *  «s  - 

and  the  horizontal  components  are  determined  by 


(11-35) 


°u  =  < 


uA  12  + 


2.3  u,  , 


mi  ' 


—  S  0 


4  >  o 


(11-36) 


The  maximum  wavelength  in  the  vertical  velocity  spectra  is 
given  by 


_ z _ 

jo  .  55  +  0.38 

5.9  z  , 


1 . 8  zi 


1  -  exp 


-4  z 


•i  ) 


( 

0 . 0003  exp 


for  z/L  <  0,  and 


0  s  zi  |L|  , 

L  <  z  <.  0 . 1  zi  , 

0.1  zi  <  z  <  zi  . 

(II-37a) 


Kv  =  z  :  Kv  *  2.9  1  ,  2  0  .  ( II~37b) 

In  the  above,  z;  denotes  the  top  of  the  planetary  boundary 
layer . 

The  parameter,  Ri,  is  the  gradient  Richardson  number, 
defined  previously  by  equation  11-19.  The  critical  Richardson 
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number,  Ric,  beyond  which  turbulence  is  suppressed,  is  given 
by 

Ric  =  0.115  Az0-175  ,  (11-38) 

where  Az  is  measured  in  centimeters.  The  Lagrangian 
autocorrelation  terms  are  determined  by 


Ru  (At) 
Rv(  At) 
i?„(At) 


(11-39) 


Tl  is  the  Lagrangian  integral  or  e-folding  time  scale  for 
velocity  autocorrelations,  and  is  determined  for  each 
component  from  the  turbulence  spectra  with 

Ttu  -  0.2P ,yv, 

Tlr-0.2  .  (11-40) 

t,.  *  c-nj-jv , 

where  Xm  is  the  dominant  wavelength  for  each  component. 
Furthermore,  the  average  velocity  is  defined  by 

V  =  )/u*  *  v2  +  w2  .  (11-41) 

the  ratio  of  Lagrangian  to  Eulerian  time  scales,  is  given 
by 
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Pv  • 


(11-42) 


P„ 


where  (3  is  restricted  to  0  <  10. 

The  horizontal  components  for  peak  wavelength  are 

given  by 


1.5zi  , 


0.7  2,, 


(11-43) 


The  constant,  A,  varies  as  a  function  of  stability  in  the 
unstable  boundary  layer,  and  is  computed  from 


A 


0.31  (1  -  3-5)  3  (1  -  15-|)  0-25  (0.55  +  0.38^)  ,  |-||sl  , 

i-t  i-t  Ll  j—l 

_  1  ~ 

0.05(1  ~  3  ~)  3(1-  15  -p )  0-25  ,  0.1  |— ij  >  \~  |  >  1  . 

■Li  i-t  Li  Li 

(11-44) 


Values  of  u'(t-At)  and  v'(t-At)  are  determined  randomly  from 
a  normal  probability  distribution  with  zero  mean.  However,  the 
turbulence  distribution  of  vertical  velocity  is  skewed  [Figure 
10],  and  the  normal  distribution  is  modified  by  changing  the 
method  of  computing  w: 


37 


aa*  +  —  -  n  ,  4  *  0  ' 

+  “  L  (11-45) 

—  +  aw  +  r|  ,  —  >  0  , 

a  L 

where 

a  =  -0.028  +  0.6 |P|  ,  (11-46) 

x\  =  0.54  \P\  , 

and 

o.i  + - — - ,  4  *  °  ' 

0.68(1  -  15-f  )  -°-25  -  1 . 8-p 

L  L 

0.1  -  0.2(  —  )0-2  ,  4  >  0  ' 

Li  Lj 

(11-47) 

The  variables  u+  and  «'  are  random  values  obtained  from  a 
normal  probability  distribution  with  zero  mean  and  standard 
deviation,  ctw.  As  noted  by  Pielke  (1984,  p.  179),  this  method 
is  not  entirely  satisfactory,  since  it  depends  on  the 
particular  random  number  generator  routine. 

3.  The  NPS  Dispersion  Model 

As  mentioned  above,  L,  the  Monin-Obukhov  length  was 
proscribed  for  this  study,  while  V,  the  mean  windspeed  at  a 
height  of  two  meters  was  set  to  3m/sec  and  the  surface 
vegetative  canopy  roughness  length  was  assumed  to  be  0.05m. 
Other  flow  parameters  reguired  by  the  McNider  and  NPS  particle 
models  were  computed  from  recently  developed  boundary  layer 
parameterizations ,  yiven  here  and  by  Panofsky  and  Dutton 
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Figure  10.  Probability  density  function  of  vertical  velocity 
in  the  convective  boundary  layer.  From  Weil,  Dispersion  in 
the  Convective  Boundary  Layer,  Lectures  on  Air  Pollution 
Modelling,  Venkatram  and  Wyngaard,  1988. 
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(1984),  Sorbjan  (1990)  and  Kamada  (1992a).  More  detailed 
descriptions  of  the  underlying  theory  will  be  published  in  a 
forthcoming  article.  Many  of  the  algorithms  listed  here  are 
actually  a  part  of  the  mesoscale  windflow  simulation  which  is 
also  used  to  drive  the  McNider  model.  From  this  windflow 
simulation  the  McNider  model  obtains  the  following:  the 
turbulent  diffusivity,  Km,  gradient  Richardson  number,  Ri, 
vertical  windshear,  the  windspeed  at  height  z,  v2,  and  surface 
layer  friction  velocity,  u.. 

u.  is  computed  from  L,  using  the  intermediate  Businger 
function,  i^m,  as  determined  by  equation  11-24,  and  the  square 
root  of  the  surface  drag  coefficient, 

,-.1/2  _  k 

^ dn  • 

so  that 


u,=MAX(C^2,  0.01)  .  (11-49) 

Given  L  and  u.,  one  can  compute  the  surface  vertical 
temperature  flux, 


= 


u.3e 

kgL 


(11-50) 


This  allows  an  estimate  of  the  free  convective  boundary  layer 
turbulent  velocity  scale  (for  L  <  0)  , 
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( 


(11-51) 


gzi 


w'ft'o 

0 


\i/3 


This  can  be  modified  to  include  shear  induced  surface  layer 
turbulence  (forced  convection)  according  to 


0.4  -  — 

1  -  m 

(l  -  150  Z°  )  -  1 

4 

l  L1/3J 

3.4 

(11-52) 


(Kamada,  1992a) . 

For  z/L  <  -0.5,  i.e.,  above  the  surface  layer,  the 
following  forms  were  used.  The  vertical  velocity  variance  was 
given  by 


o2  =  w2 


(  ^  \2nt  „  \ 

2/3 

t  ~  \2/3 

(  r*  Nl/3 

1 .  if— I  (l  -  — 

+  r2/ 3 

1  -  +  D\ 

\zi)  l  ziJ 

\zi) 

{  zi  1 

(11-53) 

where  R  =  0 . 2  and  D  =  0.1  are  ratios  of  the  inversion/surface 
temperature  flux  and  entrainment  zone/boundary  layer  depth 
(Sorbjan,  1990) . 

In  a  standard  atmospheric  windflow  model  .with  second 
order  turbulence  closure  scheme,  the  turbulence  kinetic  energy 
(tke)  is  computed  prognostically .  If  crw2  is  supplied  as  above, 
the  horizontal  variances  follow  by  subtraction.  Without  such 
a  model,  au2  and  av2  are  obtained  diagnostically  as  follows. 

The  Andre  (1978)  third  order  turbulence  closure  used 
the  following  prognostic  form  for  vertical  velocity  variance, 
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(11-54) 


3^e 

Assuming  steady  state  and  neglecting  diffusion,  this  can  be 
truncated  to  estimate  the  anisotropy,  aw2/e  where  e  is  the  tke. 
Mason  and  Thompson's  (1987)  large  eddy  simulation  of  neutrally 
stable  flow  showed  that  aw2/e  =  1/3  near  the  surface  and 
increased  with  height.  This  occurs  because  the  surface 
restricts  vertical  motion.  Also,  most  boundary  layer 
turbulence  results  from  "surface  no-slip"  induced  shear  which 
creates  mostly  horizontal  tke.  The  increase  with  height 
occurs  because  the  shear  drops  rapidly,  so  uw2/e  gradually 
approaches  the  isotropic  2/3  value  through  pressure  re¬ 
distribution,  consistent  with  Mason  and  Thompson's  results. 

Like  the  tke  on  which  it  depends,  e,  the  molecular 
dissipation  rate  of  tke  also  decays  gradually  with  height.  For 
convective  boundary  layers,  Sorbjan  (1990)  parameterizes  e  as, 

3(1  -  Z/ZJ 

0e  =  -  +  Rz/Zj  ,  and  (11-55) 

4 

e  =  06W.//Z,  .  (11-56) 

Putting  the  above  together,  the  tke  anisotropy  ratio  can  be 
parameterized  as, 
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<V  2 e0  -  e  3w'e'z 

-  =  -  +  -  ,  (11-57) 

e  3  e0  e 

(Kamada,  unpublished) .  The  first  term  on  the  right  hand  side 
accounts  for  the  height  dependence  under  neutral  stability, 
while  buoyancy  flux  in  the  second  term  accounts  for  non¬ 
neutral  stability.  This  form  actually  corresponds  to  the 
Lenschow  et  al.  (1980)  measurements  better  than  the  Andre 
model  itself  or  derivatives  thereof  (Therry  and  Lacarrere, 
1980)  . 

So  the  horizontal  turbulent  velocity  variances  become 

a2  =  (5/9)  aj(  2e/aJ  -  1)  ,  and  (II-58a) 

a2  =  O.8ou2  .  ( II-58b) 

This  allows  one  to  compute  e,  the  turbulence  kinetic  energy 
(tke)  as, 


®  =  (  <V  +  +  <V  )/2  .  (11-59) 

The  molecular  tke  dissipation  length  scale  is  parameterized  as 
1,  =  1/(1/Z  +  1.4ee1/2  )  (11-60) 

from  which  the  turbulent  mixing  length  is  estimated  as, 

lk  =  MIN  (31,  oj/e,  z)  .  (11-61) 
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This  finally  yields  the  required  turbulent  diffusivity, 


Km  =  0.44lkel/ 2  ,  (11-62) 

(Kamada,  1992b) .  The  buoyancy  gradient  for  the  unstable 
boundary  layer  is  given  by 

3dQ/dz  =  0.6Q./zi((l-z/zi)V3/  (z/zj1'3  - 

2RV3 (z/z,)2'3/  (1  -  z/zf  +  d)V3) 

(11-63) 

(Sorbjan,  1990) .  The  above  formulations  were  applied  to  the 
convective  (unstable)  boundary  layer  (z/L  <  -0.5)  above  the 
surface  layer. 

For  the  unstable  surface  layer  (L  <  0, 


but  z/L  >  -0.5), 

0m  =  (1  -  28Z/L)1'4  ,  (11-64) 

Km  =  ku.z/<pm  ,  (11-65) 

e  =  (5.6KJZJ2  ,  (11-66) 

Ri  =  z/L  ,  (11-67) 

0,  =  (  12  -  0.5ZJL  )U3  ,  (11-68) 

02  =  0,  ,  (11-69) 

03  =  (  1  -  14Z/L  ) 1/4  ,  (11-70) 

0.  =  (1  +  .  75*  (  —  Z/L) 2/3 ) 3/2  ,  (H-71) 
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and  finally 

e0  =  n.i0jzi  ,  where  the  subscript  0  refers  vlI-72) 

to  the  surface  value.  (Sorbjan,  1990) . 

For  the  neutral  to  stable  boundary  layer  and  surface 


layer,  where  L  £  0, 

a  =  2  -  10 /L  ,  and  (11-73) 

3=3-  20 /L  ,  (11-74) 

(Kamada,  1992b) .  So  that, 

0,  =  2.2(1  -  z/zt) aJ1  ,  (11-75) 

02  =  0,  ,  and  (11-76) 

0S  =  1.6(1  -  Z/Zi)a/2  .  (11-77) 


The  local  friction  velocity  and  Monin-Obukhov  lengths 
were  characterized  in  Sorbjan  (1990)  as 

u,  =  u.(l  -  z/z{) 0/2  ,  and  (11-78) 

L,  =  L  (1  -  Z/ZJ 3aJ2p  .  (11-79) 

and  the  temperature  flux  at  height,  z,  was  given  by 

w'Q'z  =  w'~Q'0(l  -  z/zjt  .  (11-80) 

The  Brunt-Vaisala  frequency  at  height,  z,  was  given  by 
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BVF  =  4 . 3u,Q.2  ( 1  +  3 . 7z /L)  (1  -  zjzj^/z 


(11-81) 


The  tke  dissipation  rate  at  z  was  parameterized  as 
<*>,  =  3.6(1  +  3 . 7z/L)  ( 1  -  z/zjt/z  ,  and 
e  =  <p(u.3 

The  Richardson  number  was  estimated  as 

z 

Ri  =  - 

5z  +  L 

and  the  momentum  diffusivity  was  given  roughly  by, 

ku.z  (1  -  z/zj 

Km  =  -  , 

1  +  5  z/L 

( Sorb j an,  1990) .  Both  the  unstable  surface  layer  and 
to  stable  boundary  layer  cases  used: 

ou  =  U.(p, 

°V  =  °u 

aw  =  u.0j 

=  2ln[  (1  +  (1  -  14z/L) in)  /2  ]  ,  and 

0;  =  9  +  9,(  ln(z/z0)  -  \ph  )  /k  ,  where 
0.  =  w'  ev9/u. 


(11-82) 

(11-83) 

(11-84) 

(11-85) 

neutral 

( II-86a) 
( II-86b) 
(II-86C) 
(11-87) 
(11-88) 
(11-89) 
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The  mean  wind  at  height  z  and  its  shear  were  estimated  to  be 
Vz  =  u.(  ln(z/z0)  -  )/k  ,  and  (11-90) 

dVz/dz  =  u.(l  +  4.7z/L)a/2(l  -  z/zt) 0/2 /kz  ,  (11-91) 

(Panofsky  and  Dutton,  1984  and  Sorbjan,  1989  ) . 

The  following  details  only  the  significant  items 
distinguishing  the  NPS  and  McNider  models.  The  NPS  particle 
model  utilized  the  above  formulations  for  the  velocity 
variances  and  tke  rather  than  those  from  McNider.  For  the  NPS 
particle  model,  the  Lagrangian  vertical  timescale  was  also 
estimated  differently  according  to 

Tlw  =  KJow  (  L  >  0  )  ,  and  (11-92) 

Tlw  =  0.3zjw.s  (  L  <  0  )  .  (11-93) 

Recognizing  that  the  horizontal/vertical  eddy  aspect 
ratio  decreases  with  increasing  stability  in  the  convective 
boundary  layer,  the  horizontal  Lagrangian  integral  timescales 
were  estimated  by 

Tlu  =  (  2  -  40 /L)  Tlw  ,  and  (11-94) 

Tlv  =  Tlu  .  (11-95) 

The  McNider  algorithms  for  horizontal  integral 
timescales  were  retained  for  stable  cases. 
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Rather  than  use  the  McNider  formulations,  the  vertical 
velocity  skewness  was  simulated  by  converting  every  fifth 
updraft  into  a  downdraft,  while  accelerating  the  updrafts  and 
decelerating  the  downdrafts.  So  the  updraft/downdraft 
probability  ratio  was  set  to  60%/40%,  while  the  updrafts  were 
50%  faster  than  the  downdrafts,  in  txme  with  atmospheric 
measurements . 

Non-stochastic  buoyant  forcing  was  also  added  to  the 
NPS  particle  model  for  non-zero  density  gradients.  For  pure 
buoyancy  driven  motion,  the  force  balance  is  just 


dw  _  _  68 

3t  9  8 


(11-96) 


where  w  is  vertical  velocity,  and  60  is  the  change  in 
potential  temperature  due  to  a  displacement,  Sz ,  away  from  the 
particle's  neutral  buoyancy  height  in  a  domain  where 
dQ/dz  0.  Here  the  neutral  buoyancy  height  is  taken  to  be 
the  initial  release  height  of  the  particle.  Then, 

dw  =  -  -f  50  at  .  (11-97) 

U 


If  the  particle1 s  "memory"  were  perfect,  then  after  a 
time,  t,  its  buoyancy  forced  velocity  would  simply  be 

t 

w  =  -  |  so  at  .  (n-98) 

o 
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However,  to  simulate  dilution  by  ambient  fluid,  the  particle 
must  gradually  forget  its  neutral  buoyancy  level.  Its 
mnemonic  e-folding  time  will  be  j?lw,  the  integral  eddy 
coherence  time  scale  already  used  to  compute  the  stochastic 
component.  So  at  each  time  step  the  particle's  memory  dims  by 
the  factor. 


M  » 


(11-99) 


In  this  case,  after  one  time  step, 

w,  =  -(g/Q)  SQjSt  ,  (Il-ioo) 

as  before.  But  for  subsequent  time  steps,  the  solutions  are 
w2  =  -(g/Q)St(MSe,  +  se2)  ,  (ii-ioi) 

w3  =  ~(g/Q)St(M2SQ1  +  M8e2  +  se3)  ,  (11-102) 

or  in  general, 

m 

=  ~(g/Q)6t  E  50,  MmJrU  .  (11-103) 

1 

This  series  grows  quickly  with  m.  However,  also  observe 
that 


K+i  =  mw,  -  (g/e)  sei+l6t  .  (II-104) 

So  only  the  last  velocity  need  be  retained.  The  buoyant 
forcing  then  added  to  the  standard  stochastic  component 
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(adjusted  for  vertical  velocity  skewness)  to  obtain  the  total 
vertical  velocity  fluctuation. 

4.  Random  Number  Generator  Kernel 

The  dispersion  models  utilize  normally  distributed 
random  numbers  to  generate  the  fluctuating  velocity.  These 
fluctuations  are  scaled  by  the  standard  deviation  of  the 
velocity  components  to  relate  them  to  the  turbulent  kinetic 
energy  of  the  flow  field.  The  method  used  to  generate  random 
numbers  in  this  study  was  a  variation  of  the  congruence 
method,  and  is  outlined  in  program  McNider,  subroutines  NORNG 
and  STRNUM  (Appendix  C)  .  To  check  the  distribution  of  the 
random  number  generation  routine,  the  count  distribution  was 
plotted  versus  standard  deviation  [Figure  11].  It  appears  to 
have  a  generally  Gaussian  shape  with  perhaps  slight  skewness 
to  the  left  of  center  for  the  3600  iterations  used  in  this 
test. 
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III.  ANALYSIS 


A.  EQUIPMENT /SOFTWARE 

The  logistics  difference  equation,  Henon  equations,  and 
the  McNider  Dispersion  model  were  all  written  in  standard 
FORTRAN-77,  and  run  on  Intel  386-based  personal  computers 
using  the  PROFORT  compiler.  All  graphs  were  produced  with 
Golden  Software's  Grapher  and  Surfer  programs. 

B.  LOGISTICS  DIFFERENCE  AND  HENON  EQUATIONS 

1.  Methodology  in  analyzing  the  Logistics  Difference 
Equation 

Computer  generated  data  from  the  logistics  difference 
equation  was  produced  using  Program  Feigenbm  (Appendix  1)  . 
Since  use  of  the  logistics  difference  equation  involves 
recursion,  an  iterative  series  rather  than  a  time  series  is 
created,  such  that  the  fractal  lengths  were  redefined  using 

T 

Hi)  =-  f  \F(i*e)  -FU)  |  di  ,  (Ill-l) 

€  J 
0 

which  is  approximated  numerically  by 

L(i)  =  ]T  | F(i  +  k)  - F(i)  |,  (III-2) 
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F(i)  =  value  of  the  function  at  the  ith  iterate, 
i  =  iteration  step, 

and  k  =  resolution  in  terms  of  number  of  iterations . 

Analogous  to  the  Lagrangian  particle  models  studied 
later,  the  resulting  iterative  series  represents  a  single 
particle,  shifting  position,  with  each  successive  iteration 
corresponding  to  a  single  time  step.  To  insure  independence 
from  initial  values,  the  program  was  run  for  3700  iterations 
and  the  first  100  points  were  discarded.  N  was  set  to  3600 
iterations  for  all  plots,  giving  3600  point  data  sets. 

Since  recursion  does  not  permit  non-integer 
iterations,  the  values  of  k  had  to  be  integer  divisors  of 
3600.  The  45  available  values  of  k  used  are  listed  as  an  input 
data  file  in  Program  Feigenbm  (Appendix  1} . 

From  a  series  of  k  values  for  L(k)  for  a  given  n ,  a 
standard  linear  regression  routine  then  determined  DA .  The 
standard  deviation  of  DA,  was  very  small  in  chaotic 
situations,  and  large  for  during  period  doubling.  The  reasons 
for  this  are  discussed  below. 

In  computing  DA,  the  last  three  values  of  L(k)  and  k 
were  discarded.  On  most  plots,  the  curve  flattened  for 
K  >  N/3.  This  was  also  noted  by  DeCaria  (1992)  in  analyzing 
time-series  data.  So  in  the  parabolic  logistics  difference 
equation,  one  reason  for  discarding  the  highest  k  values  is 
the  rising  number  of  fold-crossings  with  larger  k.  The 
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"particle"  trajectory  is  confined  to  the  unit  interval  by  the 
left-right  symmetry  or  fold  at  mid-parabola.  Therefore,  the 
distance  traversed  during  a  fold  crossing  is  not  fully 
accounted  as  with  a  system  without  folds.  The  change  in 
"apparent"  traversed  distance  will  be  lessened  as  more  fold 
crossings  are  involved.  This  results  in  decidedly  non-linear 
slopes  for  log  L(k)  versus  log  k  for  k  which  are  a  sizeable 
fraction  of  the  window. 

2.  Analysis 

For  values  of  /i  corresponding  to  a  fixed  point,  the 
plot  of  Log10L(k)  vs  Log10(k)  is  fairly  flat,  and  the 
corresponding  DA  value  is  small  [Figure  12].  For  values  of  n 
corresponding  to  period  doubling  (period  2) ,  the  plot  shows 
large  jumps  from  a  baseline,  corresponding  to  multiple 
harmonics  of  2  [Figure  13].  The  baseline  corresponds  to  a  zero 
length.  For  these  plots,  the  zero  was  adjusted  by  adding  one 
to  avoid  a  baseline  at  -  »,  a  result  of  Log10(0)  .  This  shift  by 
one  unit  was  important  for  later  calculations;  without  it,  DA 
-  *  for  all  periodic  functions. 

A  period  4  plot  still  shows  regularity  [Figure  14]. 
There  appear  to  be  three  sets  of  overlapping  slope  patterns: 
one,  a  baseline  of  slope  zero;  the  other  two  with  similar 
slopes  but  different  amplitudes.  Again,  this  can  be  ascribed 
to  harmonics.  The  length  is  zero  for  every  fourth  data  point, 
and  so  will  remain  zero  for  harmonics  having  k  values  of  4,  8, 
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Figure  12.  Logistics  Difference  Equation,  xn+I  = 
x0  =  0.01,  Fixed  Point  Attractor. 


2  •  8xn  ( l~xn)  , 


55 


Log  10  L(k) 


Lo  g  1 0  ( k ) 


Figure  13.  Logistics  Difference  Equation.  xn+1  =  3.0xn(l-xn), 
x0  =  0.01,  Period  2. 
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16,  etc.  The  length  is  nonzero  for  every  second  data  point, 
not  including  multiples  of  four.  There  is  similarity  in 
lengths  for  k  values  of  2,  6,  10,  18,  etc.  Likewise,  there  is 
similarity  in  lengths  for  k  values  which  are  not  multiples  of 
2,  e.g.,  1,  3,  5,  etc. 

At  the  onset  of  chaos  (ju  =  3.569946)  as  in  Figure  15, 
all  L(k)  lift  off  the  baseline;  without  periodicity,  all 
lengths  will  remain  non-zero.  Thus,  DA  will  increase  while  aDa 
decreases.  Also  noted  is  a  sensitivity  to  initial  conditions 
for  identical  values  of  n,  another  characteristic  of  chaos. 
Figure  16,  with  x0  =  0.04,  is  different  from  Figure  14  with  x0 
=  0.1.  However,  Figure  15,  with  x0  =  0.1,  is  virtually 
identical  to  Figure  17,  with  x0  =  0.01.  There  may  be  some 
intrinsic  pattern  to  this  sensitivity  to  initial  conditions. 
The  fold  crossing  patterns  might  provide  further  insight. 
However,  this  question  strays  from  the  present  focus. 

For  full  chaos,  aDt  becomes  relatively  small,  and  the 
plot  becomes  almost  linear,  as  shown  in  Figure  18.  This  plot 
highlights  why  the  last  three  points  were  discarded  in 
computing  DA;  the  plot  remains  very  linear  sans  large  values 
of  k.  The  main  reason  is  the  fold  crossing  patterns  mentioned 
earlier.  Additionally,  not  enough  data  points  are  retained  in 
the  length  computation  for  large  values  of  k  to  maintain 
similarity.  For  k=1200  and  N=3600,  the  computed  length  is  the 
sum  of  only  three  distance  measurements,  which  is  not  enough 
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Figure  14.  Logistics  Difference  Equation,  xn+1  =  3 . 48xn(l-xn)  , 
x0  =  0.01,  Period  4. 


58 


Log  10  L(k) 


LoglO(k) 


Figure  15.  Logistics  Difference  Equation, 

xnM  =  3 . 56994568x„(  l-xn)  ,  x0  =  0.1,  Onset  of  Chaos. 
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Figure  16.  Logistics  Difference  Equation, 

xn  +  1  =  3 . 56994  568xn  ( l-xn'  ,  x0  =  0.04,  Onset  of  Chaos 
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Figure  17.  Logistics  Difference  Equation, 

xn*,  =  3 . 56994568xn(  l-xn)  ,  x0  =  0.01,  Onset  of  Chaos. 


to  develop  a  good  statistical  mean  length. 


A  plot  of  Da  versus  n  shows  many  of  the  features  noted 
in  the  previous  bifurcation  plot  [Figure  19].  At  n  a  3.0  DA 
jumps  from  zero  (corresponding  to  a  fixed  point)  to  a  0.58, 
corresponding  to  period  2.  At  /i  «  3.45  is  seen  a  transition  to 
period  4.  Higher  values  display  regions  of  periods  8  and  16. 
However,  at  this  level  of  resolution  period  32  or  higher  order 
bifurcations  cannot  be  discerned.  Within  the  chaos  regime, 
about  ju  a  3.6,  DA  seems  to  oscillate,  and  shows  several  sharp 
peaks  and  dips.  The  apparent  window  at  n  a  3.85  corresponds  to 
a  period  4  oscillation. 

A  plot  of  entropy  versus  n  shows  similar  features 
[Figure  20].  Regions  of  period  2,  4,  and  8  are  readily 
apparent.  As  in  the  DA  versus  n  plot  [Figure  19],  differences 
between  period  16  and  above  and  chaos  cannot  be  discerned  at 
this  level  of  resolution.  The  large  window  of  low  periodicity 
at  /x  a  3.85  is  also  seen,  as  are  several  other  windows  of 
periodicity  at  roughly  n  a  3.63  and  n  a  3.73. 

The  Lyapunov  exponent  also  shows  the  trends  seen  in 
the  fractal  dimension  and  entropy  [Figure  21].  X  dips  sharply 
at  4  a  3.25,  3.5,  and  3.55,  the  centers  for  periods  2,  4,  and 
8.  Again,  changes  beyond  period  16  are  hard  to  resolve.  In  the 
period  doubling  region,  X  =  0  indicates  that  the  function  is 
transitioning  to  the  next  bifurcation,  a  fixed  point 
fissioning  to  two  fixed  points.  For  X  >  0,  the  function  is 
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Figure  18.  Logistics  Difference  Equation, 
xn  +  1  =  3 . 999 6xn ( l-xn)  ,  Xq  =  0.01,  Chaos. 
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Figure  20.  Shannon  Entropy  Analysis  of  Logistics  Difference 
Equation,  xn+1  =  MXn(l-xn). 
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Figure  21.  Lyapunov  Exponent  (X)  Analysis  of 
Difference  Equation,  xn+i  =  jixn{l-xJ. 


Logistics 
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chaotic.  There  are  also  numerous  "windows  of  periodicity"  in 
the  chaotic  region  which  appear  in  the  bifurcation,  fractal 
dimension,  and  entropy  plots.  These  windows  are  more  readily 
identified  by  the  Lyapunov  exponent;  if  X  drops  below  zero, 
the  function  is  periodic. 

Figure  22  shows  the  strong  correspondences  between 
entropy,  fractal  dimension,  and  Lyapunov  exponent. 
Periodicities  and  chaos  are  apparent  with  each  metric,  as  are 
some  differences.  DA  sometimes  peaks  where  S  and  X  dip,  as  at 
M  ~  3.725.  The  low  values  of  X  and  S  indicate  periodic 
regions,  for  which  DA  sometimes  sharply  increases.  The 
bifurcation  map  [Figure  6]  provides  no  clue  as  to  the  cause. 
However,  if  the  fixed  points  are  more  widely  separated  than 
the  average  distance  between  successive  points  in  the 
surrounding  chaos  regions,  DA  will  become  relatively  large. 
This  is  less  likely  for  real  fluids  because  motion  is  then 
constrained  by  physical  forces  and  conservation  laws.  However, 
it  points  to  a  possible  schism  between  diffusion  and 
dispersion  rates  even  without  mean  flow. 

3.  Henon  Function  Analysis 

Since  the  Henon  equations  are  a  2-D  extension  of  the 
logistics  difference  equation,  similar  behavior  is  expected. 
With  two  dimensions,  the  length  measurement  is  modified  to 
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Figure  22.  Fractal  Dimension  (DA)  ,  Shannon  Entropy  (S) ,  and 
Lyapunov  Exponent  (X)  for  Logistics  Difference  Equation. 
Left  axis  X;  right  axis  S  and  DA. 
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Ar  =  sjU^-xJ2  +  (ya#1  ~yn) 2 


(III-3) 


In  the  period  2  case  [Figure  7],  a  plot  of  LogI0L(k)  vs 
Log10k  is  very  similar  to  the  logistics  difference  equation  for 
period  2  [Figure  23].  Again,  for  every  k,  evenly  divisible  by 
2,  the  length  drops  to  zero.  There  are  two  sets  of  slope 
patterns:  a  baseline  with  slope  zero,  and  an  upper  non-zero 
slope  (corresponding  to  non-even  k-values) .  This  upper  slope 
appears  to  be  quite  straight,  except  for  the  largest  k  value. 

At  the  onset  of  chaos  [Figure  24],  again  the  baseline 
lifts,  and  the  overall  deviation  in  slope  is  smaller 
[Figure  25].  It  is  also  noted  that  Log10L(k)  only  approaches 
zero  when  k  is  quite  large.  As  noted  before,  for  large  k  there 
are  not  enough  distance  measurements  to  adequately  define  a 
good  statistical  mean  length;  nor  are  there  enough  points  to 
adequately  define  a  periodicity. 

At  full  chaos  [Figure  7],  the  plot  of  Log10L(k)  vs 
Log10k  is  so  linear  that,  in  regions  of  full  turbulence  or 
chaos,  Da  can  be  defined  with  only  two  date  points,  i  and  o, 
corresponding  to  one  small  inner  and  one  large  outer  scale  of 
k.  [Figure  26].  Then  DA  in  this  case  is  given  by 
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Figure  25.  Henon  Function,  a  =  1.057,  b  =  0.3,  DA  =  1.114, 
Onset  of  Chaos. 
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This  shows  immediately  that  an  increase  in  DA  implies  a 
relative  increase  in  the  apparent  magnitude  of  the 
fluctuations  discernible  at  higher  resolutions,  i.e.,  a  shift 
in  the  amplitude  spectrum  toward  smaller  scales.  This 
interpretation  of  DA  will  be  of  use  in  analyzing  the 
atmospheric  Lagrangian  particle  models. 

There  appears  to  be  a  close  correspondence  between 
entropy  S  and  DA  for  the  Henon  function  [Figure  27].  Both  show 
jumps  corresponding  to  periods  2,  4,  8,  and  perhaps  even  for 
period  16.  Both  plots  show  corresponding  changes  when  there  is 
periodicity  within  the  chaos  region. 

A  higher  resolution  look  at  the  region, 
1.052  <  a  <  1.082,  displays  some  interesting  complexity 
[Figure  28].  Not  only  are  there  normal  bifurcations,  but  at 
a  as  1.062  a  new  branch  appears  with  no  obvious  connection  to 
previous  points.  This  does  not  correspond  to  the  fissioning 
process  observed  earlier,  but  rather  to  an  entirely  new  set  of 
solutions.  At  a  »  1.080,  points  from  this  new  branch  break 
away  and  drift  toward  the  original  branch. 
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In  the  region,  1.052  <  a  <  1.082,  DA  jumps  in 
amplitude  whenever  the  bifurcation  map  shows  a  window  of 
periodicity  [Figure  29].  DA  based  on  x  rather  than  the  total 
distance  r  portrays  much  but  not  all  the  same  information: 
jumps  correspond  to  most  of  the  same  windows  of  periodicity 
[Figure  30].  The  strong  peaks  in  DA  correspond  to  sharp  dips 
in  S;  however,  sharp  dips  in  S  do  not  necessarily  correspond 
to  sharp  peaks  in  DA  [Figure  31].  So  a  jump  in  average 
distance  traversed  between  iterations  always  accompanies  a 
drop  in  position  randomness  but  not  vice  versa.  This  suggests 
that  periodicity  in  the  Henon  system  sometimes  results  from 
the  sudden  appearance  of  fixed  points  which  are  as  closely 
spaced  or  more  closely  spaced  than  the  mean  distance  between 
successive  iterations  in  the  surrounding  chaos. 

Recognizing  that  there  are  two  free  parameters  in  the 
Henon  equation  (a  and  b) ,  DA  and  S  were  mapped  for  values  of 
0  <  a  <  1.5  and  0  <  b  <  1.0  [Figures  32,  33].  The  two  3-d  maps 
are  remarkably  similar,  and  suggest  that  exp(DA)  would  be  of 
the  same  order  as  S.  The  regimes  of  low  periodicity  are  well 
defined.  Where  DA  and  S  are  zero,  the  Henon  function  has 
period  1,  corresponding  to  one  fixed  attractor.  The  first 
jumps  to  periods  2  and  4  are  well  defined  at  this  level  of 
resolution.  There  are  other  steps  of  discernible  width  for 
high  a  and  b  values,  on  the  left  and  back  sides  of  the 
graphical  "mountain".  Another  interesting  feature  is  the 
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Figure  29.  Fractal  Dimension  ( DA)  Analysis  of  Henon 
Function,  testing  Ar  with  b  =  0.3. 
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Figure  30.  Fractal  Dimension  (DA) 
Function,  testing  Ax,  with  b  =  0.3. 


behavior  for  high  values  of  b,  where  both  DA  and  S  show  sudden 
increases  in  value. 

4.  Conclusions  drawn  from  Logistics  Difference  and  Henon 

equation  analysis 

The  self-affine  fractal  dimension,  DA,  is  a  good 
discriminator  of  low  frequencies  in  data,  e.g.,  periods  2,  4, 
and  8  of  a  given  length  data  set.  However,  at  high 
frequencies,  period  16  and  above  in  the  initial  bifurcation, 
it  is  difficult  to  use  DA  to  differentiate  between  periodicity 
and  turbulence  or  chaos.  The  entropy,  S,  seems  related  to  DA, 
but  S  measures  the  evenness  of  particle  position  distribution, 
while  Da  measures  the  changing  apparent  jaggedness  of  the 
function  as  the  resolution  varies.  A  sharp  change  in  S  always 
accompanies  a  sharp  change  in  DA ,  but  not  always  vice  versa. 
Therefore,  DA  and  S  generally  show  the  same  trends  within  the 
chaos  region,  but  not  always.  In  such  cases  diffusion  and 
dispersion  rates  are  not  equivalent. 

The  Lyapunov  exponent  provides  perhaps  the  most  definitive 
information:  for  a  given  value  of  X,  we  know  for  certain 
whether  the  function  is  stable,  periodic,  or  chaotic. 
Additionally,  the  Lyapunov  exponent  should  be  able  to  describe 
the  degree  of  chaos:  the  greater  the  value  of  X,  the  faster 
the  orbital  trajectory  diverges,  and  therefore  the  more 
chaotic  the  function.  Unfortunately,  the  Lyapunov  exponent  is 
not  readily  determined  for  multi-dimensional  systems. 
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Figure  32.  Fractal  Dimension  (DA)  Analysis  of  Henon 
Function,  xn  +  1  =  l-axn2+yn,  yn+1  =  bxn. 
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Figure  33.  Entropy  (S)  Analysis  of  Henon  Function, 
=  i-axn2-yn»  yn+.  =  bxn. 
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Analysis  of  these  two  simple  systems,  the  logistics 
difference  and  Henon  sets,  shows  that  S,  DA ,  and  X  display 
correlated  but  not  identical  behavior,  since  they  measure 
rather  different  things.  These  insights  can  now  be  applied  in 
studying  atmospheric  Lagrangian  particle  models. 

C.  MCNIDER  PARTICLE  DISPERSION  MODEL  ANALYSIS 
l.  Methodology 

Four  aspects  of  the  McNider  Particle  Dispersion  Model 
were  studied:  the  divergences  of  vertical  position,  total 
position,  velocity,  and  phase  velocity.  The  position  was 
measured  in  radial  distance  from  the  particle  origin,  which 
was  arbitrarily  set  at  x,y  =  0,  z  =  100  meters.  To  simulate 
the  real  atmosphere,  the  boundary  layer  inversion  height  was 
varied  linearly  from  2,000  to  200  meters  as  1/L,  the  inverse 
Monin-Obukhov  length,  was  varied  from  -0.2  to  0.1.  Mean 
velocities  were  set  to  zero  and  total  particle  velocity 
reflection  was  assumed  at  the  ground  surface  and  inversion. 

Hints  as  to  model  behavior  are  again  provided  by  the 
bifurcation  maps,  which  now  depict  the  expansion  and 
distribution  of  particle  range  with  decreasing  stability,  as 
measured  by  1/L.  Model  performance  was  checked  against 
standard  geophysical  measures  such  as  the  Brunt  Vaisala 
Frequency  ( BVF) ,  turbulent  kinetic  energy  (tke) ,  and  vertical 
velocity  variance  (aj)  ,  as  well  as  against  the  two  readily 
calculated  chaos  measures,  entropy,  S,  and  the  self-affine 
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fractal  dimension,  DA,  for  the  range  of  1/L  values  seen  in  the 
atmosphere.  Using  100  partitions,  the  computed  maximum 
possible  entropy  for  all  plots  in  the  analysis  was 

Sm «  -4.6  .  (III-5) 

For  this  study  1/L  was  set  at  a  particular  value,  a 
single  particle  was  released  and  followed  for  3,600  time 
increments  of  1/6  seconds  each  (600  seconds  total) ,  then  1/L 
was  incremented  and  the  analysis  repeated,  until  the  entire 
range  was  spanned  using  1/L  increments  of  0.01. 

2.  McNider  No-Skew  Routine 

For  real  unstable  atmospheres,  the  vertical  velocity 
distribution  is  negatively  skewed  [Figure  10] .  Updrafts  have 
higher  velocities  and  thus  occupy  less  volume  than  downdrafts. 
However,  for  the  initial  analysis  of  the  McNider  model,  the 
vertical  velocity  turbulence  distribution  was  not  skewed.  This 
was  done  to  check  the  model  without  the  ad  hoc  method  McNider 
developed  to  introduce  skewness,  and  which  as  outlined  in 
Pielke  (1984,  pp.  178-179)  appears  to  be  incorrect. 

An  r  versus  1/L  plot  of  the  particle  position  shows 
that  on  the  negative  (unstable)  side,  DA  and  S  follow  roughly 
the  same  trend  up  to  1/L  =  0  [Figure  34].  From  the  analysis  of 
Da  in  the  Henon  system,  this  suggests  that  the  ratio  of  large 
to  small-scale  vertical  distances  between  points  varies  little 
in  the  unstable  regime  without  overlaying  skewness.  However, 
Da  again  jumps  suddenly  when  1/L  exceeds  zero.  Since  vertical 
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fluctuations  tend  to  decay  in  general  with  increasing 
stability,  one  expects  that  S  would  decay  and  DA  would  tend 
upward  steadily  from  left  to  right  for  real  turbulence  in  the 
atmosphere.  Again,  the  discontinuities  near  1/L  =  0  are 
probably  due  to  discontinuities  across  the  neutral  transition 
in  the  McNider  algorithms.  The  small  scale  fluctuations  in  S 
with  1/L  are  probably  due  to  the  random  nature  of  the  particle 
paths . 

On  the  positive  (stable)  side  of  1/L,  both  DA  and  S 
show  sudden  dramatic  step  increases,  followed  by  DA  tending 
upward  while  S  tends  downward.  Again,  this  highlights  the  fact 
that  S  and  DA  are  not  measuring  the  same  thing.  The  sudden 
jumps  in  S  in  the  McNider  model  seem  antithetical  to  the  fact 
that  negative  buoyancy  tends  to  suppress  vertical  motion  and 
render  it  oscillatory  in  real  fluids.  For  growing  positive 
values  of  1/L,  the  decrease  in  S  shows  that  the  distribution 
of  points  within  the  domain  becomes  less  uniform,  while  the 
rise  in  DA  suggests  that  the  distance  between  successive 
points  decreases  more  slowly  on  the  small  scale  than  for  the 
larger,  longer  time  scale  fluctuations.  This  is  qualitatively 
consistent  with  measurements  under  increasingly  stable 
conditions  in  real  fluids  (Stull,  1988) . 

The  bifurcation  diagram  [Figure  35]  also  shows  that 
for  1/L  <  0  the  total  range  of  positions  steadily  decreases 
with  increasing  stability,  while  for  1/L  >  0  it  hardly  varies. 
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This  is  because  the  total  position  change  includes  both 
vertical  and  horizontal  fluctuations,  and  horizontal 
fluctuations  are  much  less  dependent  on  stability. 

These  plots  reveal  other  weaknesses  in  the  McNider 
model.  Unrealistically  sharp  fluctuations  seem  to  occur  near 
the  neutral  stability  transition  (1/L  =  0),  probably  due  to 
the  dichotomous  nature  of  the  McNider  algorithms,  one  set  for 
stable  conditions,  another  for  unstable.  From  equation  11-18, 
it  is  also  clear  for  isotropic  turbulence  that  aw2  =  2/3  tke. 
Stability  suppresses  vertical  turbulence,  while  instability 
enhances  it.  So  for  1/L  <  0,  aw2  should  exceed  2/3  tke.  Figure 
34  shows  increasingly  qualitatively  incorrect  ratios  with 
increasing  instability,  aj  also  drops  to  near  zero  at  neutral 
stability  (1/L  s  0)  due  to  very  low  mean  windshear  values  in 
the  mesoscale  flow  simulation.  This  may  not  be  as  significant 
a  problem  in  the  prognostic  windflow  models  for  which  the 
McNider  model  was  originally  designed  as  it  is  in  the 
similarity-based  flow  model  simulator  used  in  this  study. 
However,  an  artificial  "dip"  would  probably  still  appear. 

The  bifurcation  map  of  the  total  3-D  position  seems  to 
show  a  fairly  constant  range  for  1/L  >  0,  a  jump  at  zero,  and 
a  varying  range  for  1/L  <  0  [Figure  35] .  The  plot  density  is 
too  heavy  to  discern  actual  distribution  patterns  for  given 
values  of  1/L,  although  the  computer  screen  displays  a  Moire¬ 
like  pattern.  What  can  be  extracted  from  this  plot  is  the 
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Figure  34.  McNider  Particle  Dispersion  Model.  Test  of  total 
3-D  distance,  Ar,  with  no  skew. 
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general  density  of  points  and  the  upper  and  lower  boundaries 
(in  magnitude)  of  position. 

In  contrast  to  Figure  34,  a  plot  of  only  the  vertical 
position  change  (Az)  reveals  that  the  pattern  for  DA  is  quite 
similar  to  that  for  S  [Figure  36].  The  bifurcation  map  shows 
an  opposite  jump  in  the  vertical  position  as  1/L  crosses  the 
zero-point  [Figure  37]. 

The  velocity  distribution  for  various  values  of  1/L 
seems  fairly  constant,  so  that  the  S  curve  has  only  a  slight 
upward  trend  [Figure  38],  except  at  values  near  1/L  =  0. 
However,  DA  tends  downward  strongly  from  left  to  right,  jumps 
at  zero,  then  again  tends  downward  at  around  1/L  =  0.7,  where 
it  begins  increasing  in  value.  This  upward  trend  in  DA  is 
probably  caused  by  increasing  negative  buoyancy.  The 
bifurcation  map  [Figure  39]  shows  that  the  velocity  range 
tends  generally  downward,  with  a  jump  near  1/L  =  0,  and  a 
change  in  trend  at  1/L  =  0.7. 

The  phase  velocity  plots  are  almost  identical  to  the 
total  position  plots.  This  is  because  in  calculation  phase 
velocity  the  velocities  are  small  in  magnitude  compared  to  the 
position  values. 

3.  McNider  Model,  Skew  On 

A  glance  at  Figure  42  shows  a  big  problem  with  the 
vertical  velocity  skewness  algorithm.  The  tke  and  ctw2  values 
are  much  too  large.  An  expected  maximum  aj  will  be  -  10m2/s2, 
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Figure  37.  McNider  Particle  Dispersion  Model  Bifurcation 
Map,  Test  of  vertical  distance,  Az ,  with  no  skew. 
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Figure  38.  McNider  Particle  Dispersion  Model,  Test  of  total 
velocity,  no  skew. 
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Figure  40.  McNider  Particle  Dispersion  Model.  Test  of  Phase 
Velocity,  with  no  skew. 
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Phase  Velocity  Range 


while  the  plot  shows  -  100  m2/s2.  It  seems  that  the  algorithm 
is  incorrect  for  very  unstable  values  of  1/L.  However,  the 
stable  side  (positive  values  of  1/L)  shows  well  behaved, 
reasonable  values  for  tke  and  crw2. 

Examining  only  the  stable  side  (positive  1/L)  of 
Figures  42-49,  DA  and  S  are  surprisingly  flat  for  velocity, 
but  not  for  total  3-D  position,  vertical  position  change,  or 
phase  velocity.  Also,  the  bifurcation  plot  of  velocity  range 
becomes  quite  small  for  positive  1/L  (Figure  47)  .  Oddly 
enough,  the  small  values  begin  not  at  1/L  >  0,  but  1/L  >  0.02. 
This  corresponds  to  a  sharp  drop  in  entropy  and  fractal 
dimension  near  this  value  of  1/L,  suggesting  that  the  skewness 
algorithm  is  responsible.  However,  re-examining  Figures  36  and 
38  with  "no-skew",  entropy  and  DA  drop  sharply  at  1/L  =  0.02 
as  well  as  at  1/L  =  0,  suggesting  that  something  other  than 
the  skewness  algorithm  is  responsible,  perhaps  the  inherent 
dichotomy  in  equations  11-44  or  II-37a. 

Again,  this  sharp  transition  near  neutral  values  of 
1/L  does  not  seem  to  accurately  reflect  nature,  but  rather 
appears  to  be  a  discontinuity  in  the  algorithm  solutions. 

Figure  47  also  shows  that  the  velocity  range  is  much 
too  high  for  the  unstable  case;  an  expected  velocity  value  is 
on  the  order  of  10  m/s  or  less.  Oddly,  there  are  several 
windows  where  the  velocity  range  drops  to  low  values,  and  are 
roughly  what  is  expected,  i.e.,  for  1/L  «  -0.475.  However,  at 
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Figure  42.  McNider  Particle  Dispersion  Model.  Test  of  total 
position,  Ar,  with  skewness. 
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Figure  45.  McNider  Particle  Dispersion  Model  Bifurcation 
Map.  Test  of  vertical  position,  Az ,  with  skewness. 
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Figure  46.  McNider  Particle  Dispersion  Model.  Test  of  total 
velocity,  with  skewness. 
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Figure  47.  McNider  Particle  Dispersion  Model  Bifurcation 
Map.  Test  of  total  velocity,  with  skewness. 
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more  negative  values  of  1/L  the  aw2/tke  ratio  indicates  that 
nearly  all  the  turbulence  is  in  the  vertical  component  which 
is  not  seen  in  the  atmosphere  nor  is  physically  likely  (Stull, 
1988)  . 

4.  NPS  Model,  Skew  On 

The  NPS  model  employed  the  double  Gaussian  skewness 
scheme  as  described  in  the  theory  section.  It  avoids  the 
problems  seen  with  the  McNider  model  skewness  scheme  regarding 
the  inappropriate  high  values  for  velocity,  tke,  and  crw2  for 
unstable  values  of  1/L.  Also,  the  tke  and  aw2  better  reflect 
the  behavior  of  real  fluids:  aw2/tke  exceeds  2/3  for  unstable 
values  of  1/L,  but  not  for  stable  values.  The  sharp 
discontinuity  in  transition  from  unstable  to  stable  1/L  is 
reduced  somewhat. 

Curiously,  there  is  a  distinct  pattern  to  the  strong 
fluctuations  in  DA,  S,  tke,  and  aw2  for  unstable  1/L  but  not 
for  stable  1/L  values.  This  was  also  observed  in  the  McNider 
skewness  routine.  The  cause  has  not  been  determined;  it  could 
be  due  to  limitations  in  the  random  number  generator  or 
perhaps  an  undiscovered  program  error.  This  is  most  obvious  in 
the  bifurcation  map  of  vertical  position  range,  Az  [Figure 
53],  where  there  are  sharp  jumps. 

The  phase  velocity  plots  are  similar  to  the  total 
position  plots  in  all  cases  examined.  The  reason  is  readily 
determined:  the  phase  velocity  is  dominated  by  the  total 
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position,  which  varies  from  near  zero  to  1000  meters,  while 
the  velocity  varies  from  near  zero  to  roughly  3  meters/ second . 

As  in  the  McNider  plots  there  is  little  indication 
that  the  entropy  is  affected  by  a  rising  BVF,  while  DA  seems 
to  rise  as  BVF  rises,  most  notably  in  the  vertical  position 
plots.  The  BVFs  indicate  a  long  period  compared  with  the 
1/6  second  timesteps.  The  entropy  analysis  will  not  detect  a 
particle  riding  on  a  low-frequency  wave  with  long  period.  In 
such  a  case  the  particle  distribution  over  time  may  appear  to 
be  roughly  uniform  between  a  minimum  and  maximum  value.  This 
constraint  is  similar  to  the  resolution  requirements  for 
entropy  computation  noted  earlier.  If  the  At  in  the  model  were 
increased  from  1/6  second  to  say  10  seconds,  while  still 
retaining  3600  iterations,  the  entropy  might  also  show  a 
variance  with  the  BVF. 
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Figure  52.  NPS  Particle  Dispersion  Model.  Test  of  vertical 
position,  Az ,  with  skewness. 
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Figure  53.  NPS  Particle  Dispersion  Model  Dispersion  Map. 
Test  of  vertical  position,  Az,  with  skewness. 


Ill 


NPS  Particle  Dispersion  Metrics 


Figure  54.  NPS  Particle  Dispersion  Model.  Test  of  velocity, 
with  skewness. 


Figure  55.  NPS  Particle  Dispersion  Model  Bifurcation  Map 
Test  of  velocity,  with  skewness. 
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Figure  56.  NPS  Particle  Dispersion  Model.  Test 
velocity,  with  skewness. 
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IV.  CONCLUSIONS 


Chaos  metrics  such  as  the  self-affine  fractal  dimension, 
Da ,  entropy,  S,  and  Lyapunov  exponent,  X,  are  useful  tools  in 
the  analysis  of  system  characteristics  and  behavior.  DA 
measures  the  jaggedness  of  a  trajectory,  or  relative  magnitude 
of  small  versus  large  scale  fluctuations.  The  entropy  measures 
the  randomness  of  state  distribution.  The  Lyapunov  exponent, 
difficult  to  implement  in  higher  dimensions,  measures  the 
divergence  of  trajectories,  and  leads  to  a  predictability  time 
scale.  As  in  both  the  chaos  and  atmospheric  systems,  DA  and  S 
often  correspond  closely,  but  not  always.  The  Henon  equations 
also  demonstrated  that  exp(DA)  is  often  of  the  same  order  as 
S. 

These  chaos  metrics  when  used  in  conjunction  with  standard 
geophysical  measures  such  as  turbulent  kinetic  energy,  tke, 
vertical  velocity  variance,  aw2,  and  the  Brunt-Vaisala 
Frequency,  BVF,  can  be  applied  usefully  to  atmospheric  3-D 
particle  dispersion  models.  The  McNider  model  as  listed  in 
Pielke  (1984)  has  an  inherent  weakness  in  the  skewness  of 
vertical  velocity  variance.  This  weakness  leads  to  obviously 
incorrect  values  of  tke,  crw2,  positions,  and  velocities  for 
unstable  values  of  the  inverse  Monin-Obukhov  length,  1/L.  The 
McNider  algorithms  also  display  an  inherent  discontinuity  as 
the  inverse  Monin-Obukhov  length  crosses  the  zero  point, 
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reflecting  a  transition  from  unstable  to  stable  buoyancy.  This 
sharp  discontinuity  does  not  correspond  to  the  real 
atmosphere,  which  suggests  that  model  performance  suffers  at 
small  values  of  1/L. 

The  NPS  particle  dispersion  model's  measured  performance 
seems  to  better  reflect  reality,  it  does  not  share  the  McNider 
model's  weakness  in  the  skewness  of  vertical  velocity 
variance,  and  displays  more  realistic  ratios  of  aw2  to  tke.  The 
NPS  Particle  Dispersion  Model  also  reduces  the  discontinuity 
in  transition  from  negative  to  positive  values  of  1/L. 
However,  the  NPS  model  still  has  a  problem  with  unstable 
inverse  Monin-Obukhov  lengths  (1/L  <  0)  ;  it  displays  a  pattern 
to  the  vertical  position  distribution  not  reflective  of  real 
fluid  behavior.  Since  the  NPS  model  was  developed  in  response 
to  the  current  study,  there  may  still  be  problems  in  coding  or 
in  the  random  number  generator.  The  NPS  model  requires  further 
refinement  to  model  particle  behavior  in  a  fluid 
realistically. 

Previous  performance  metrics  for  atmospheric  particle 
models  were  only  designed  to  measure  aggregate  particle 
behavior.  The  above  chaos  metrics  parameters  also  offer  some 
insight  into  the  model  behavior  of  individual  particles. 
Results  indicate  that  atmospheric  dispersion  models  require 
further  development  at  a  fundamental  level  in  replicating 
turbulent  diffusion.  For  example,  large  fractal  dimensions  in 
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atmospheric  Lagrangian  particle  models  seem  indicative  of 
strongly  stable,  rather  than  unstable,  conditions,  contrary  to 
what  seems  likely  for  real  turbulence.  This  points  to  the 
neglect  in  such  models  of  motions  at  scales  other  than  the 
dominant  scale  determined  by  auvw.  Thus,  atmospheric 
Lagrangian  particle  models  may  simulate  only  single  scale  or 
at  most  a  limited  range  of  large  scale  diffusion  processes. 
Fluctuations  more  in  accord  with  real  velocity  spectra  might 
display  more  realistic  entropy  and  fractal  behavior.  As  is,  S 
and  Da  behave  oppositely  at  times,  which  indicates  that 
diffusion  and  dispersion  rates  do  not  have  equivalent  meaning, 
even  in  the  absence  of  mean  windflow. 

Simple  periodic  behavior  preliminary  to  chaos  in  classical 
bifurcatory  systems  is  also  apparently  not  entirely  equivalent 
to  laminar  wave  behavior  prior  to  the  onset  of  turbulence.  For 
example,  if  the  time  resolution  of  the  analysis  is  much 
shorter  than  the  period,  the  Shannon  entropy  for  an 
atmospheric  Laa^angian  particle  model  may  be  quite  large  for 
wave  motion,  but  small  for  periodic  behavior.  This  is  because 
the  number  of  possible  states  in  wave  motion  is  not  limited  to 
the  amplitude  extrema  as  in  the  periodic  motion. 

One  weakness  of  this  study  of  Lagrangian  particle  model 
performance  is  the  unavailability  of  real  data  for  comparison 
purposes.  As  stated  earlier,  atmospheric  data  is  normally 
obtained  in  the  Eulerian  rather  than  Lagrangian  frame.  This 


118 


suggests  a  need  for  development  of  Lagrangian-based  sensing  in 
fluid  turbulence  experiments.  Position  and  velocity 
measurements  of  particle  markers  in  a  wide  range  of 
atmospheric  conditions,  both  stable  and  unstable,  as  a 
function  of  time,  with  rapid  response  remote  sensors,  would  be 
useful  in  developing  better  models  of  particle  dispersion  and 
diffusion.  As  yet,  since  the  gathering  of  data  from  real 
fluids  in  the  Lagrangian  frame  is  not  feasible,  these  study 
results  suggest  that  more  extended  application  of  chaos 
metrics  to  Eulerian  measurements  of  turbulence  in  real  fluids 
is  warranted  in  order  to  gauge  the  performance  of  Eulerian 
turbulence  models.  Such  applications  need  not  be  restricted  to 
small-scale  phenomena,  since  mesoscale  Rayleigh-Benard 
convection  also  displays  transitions  from  laminar  cellular  to 
chaotic  behavior  (Agee  et  al . ,  1973).  Model  improvements  based 
on  such  tests  may  then  be  extended  to  the  Lagrangian  frame. 

With  regard  to  the  Eulerian  frame,  both  fractal  dimension 
and  Shannon  entropy  are  simple  calculations  which  can  be 
performed  real-time  in  situ  with  current  sampling  devices,  and 
should  be  simpler  to  implement  than  FFTs ,  since  there  are  no 
transform  matrices.  One  area  to  consider  when  conducting 
measurements  of  the  self-affine  fractal  dimension  is  the  size 
of  both  e,  the  time  increment,  and  T,  the  time  window  width 
over  which  the  length  is  being  defined,  e  has  a  minimum  size 
dictated  by  the  response  time  of  the  sensors,  while  the  size 
of  T  is  critical  when  looking  for  intermittent  or  low 
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frequency  events,  e.g.,  gravity  waves  governed  by  the  Brunt- 
Vaisala  Frequency.  If  emin  is  too  small,  then  the  time 
increment  will  fall  in  the  range  where  the  velocity 
autocorrelation  decay  is  not  exponential  but  constrained  to  be 
parabolic  due  to  continuity  and  the  viscosity  of  real  fluids. 
At  the  same  time  e,,^  must  be  at  least  -  3  orders  of  magnitude 
smaller  than  T  to  establish  statistical  validity.  Establishing 
the  appropriate  is  crucial  to  distinguish  waves  from 
turbulence,  and  probably  applies  equally  for  distinguishing 
intermittent  and  coherent  structures  in  general.  (Kamada  and 
DeCaria,  1992) 

Another  area  for  further  study  is  using  the  Lyapunov 
exponent  in  analyzing  3-D  turbulence.  Although  difficult  to 
compute  in  3-D,  X  provides  a  definite  test  of  chaos  in  the 
particle  diffusion  rate.  A  study  of  all  three  chaos  metrics 
seems  necessary  to  clarify  the  relationship  of  diffusion  rate 
to  dispersion  rate  in  these  models. 

Standard  geophysical  turbulence  measures  such  as  tke  and 
aj  might  also  be  applied  to  classical  chaotic  systems  such  as 
the  logistics  difference  and  Henon  equations.  The  definitions 
of  timescale  and  averaging  time  must  first  be  established, 
e.g.,  one  iteration  equals  one  time  unit.  The  corresponding 
tke  would  be  zero  for  a  fixed  point,  and  definite  magnitude 
changes  would  occur  as  the  function  bifurcates  to  periods  2, 
4,  8,  and  to  chaos. 
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The  utility  of  chaos  metrics  in  analyzing  the  3-D 
Monte  Carlo  based  particle  dispersion  models  also  suggests 
possible  utility  for  other  laminar/turbulent  phenomena.  These 
might  include  plasmas,  acoustics,  and  laser  cavities.  Solid 
state  free  electron  gas  systems  may  also  find  these  chaos 
metrics  useful  in  describing  phonon  and  electron  transport. 
Phonon  resonance  in  crystal  lattices  also  reminds  one  that 
particle  behavior  in  lattice  gases  and  cellular  automata  may 
be  studied  with  such  chaos  metrics.  These  measures  might  also 
be  useful  in  modeling  gas  or  liquid  phase  complex  chemical 
kinetics  or  radiation,  which  have  long  been  simulated  by  Monte 
Carlo  schemes. 
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APPENDIX  A.  LOGISTICS  DIFFERENCE  EQUATION  ANALYSIS  PROGRAM 


PROGRAM  Feigenbm 
c 

c  This  program  computes  self-affine  fractal  dimension,  Da,  for  the 
c  Feigenbaum  recursion  relation, 
c 

c  x(i+l)  =  4  lambda  x(i)(  1  -  x(i), 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


as  a  function  of  lambda.  Da  is  defined  for  the  iteration  series  as: 


Da  =  - 


d  Ln[  L(i)  ] 
d  Ln[  i  1 


f  T 


,  where  L(i)  =  (1/i) 


|  F ( i+1 )  -  F ( i )  di 


N 

=  SUM  | F ( i )  -  F(i-l) 
i=k,  k 


and  k  is  an  exact  integer  divisor  of  N.  Here  we  choose 
N  =  3,600  to  give  us  at  least  three  decades  to  use  in  computing  Da. 
This  also  allows  48  exact  integer  divisors  of  N  in  the  k  array.  So 
this  gives  us  48  points  of  i  vs.  L(i)  in  the  linear  regression  for  Da. 


DOUBLE  PRECISION  lambda,  lglngth(48,  100),  lgk(48),  y,  sum,  x 
DIMENSION  y ( 0 : 3600 ) ,  sum(48,  100),  k(48),  nyval(48),  xy(2,48) 

&  ans(16),  Dalmbda(3, 100) 

DATA  k/1,  2,  3,  4,  5,  6,  8,  9,  10,  12,  15,  16,  18,  20,  24,  25, 

&  30,  36,  40,  45,  50,  60,  72,  75,  80,  90,  100,  103,  106,  109, 

&  116,  120,  144,  150,  180,  200,  225,  240,  300,  360,  400,  450, 

&  600,  720,  900,  1200,  1800,  3600/ 

OPEN  (1,  FILE=' feigenis.dat' ,  RECL=255,  STATUS= *  new ’ ,  ERR=700) 
OPEN  (2,  FILE=' feigenll.dat' ,  RECL=255,  STATUS= ' new’ ,  ERR=800) 
c 

y ( 0 )  =0.0 

c  Get  Da  for  different  values  of  lambda  from  0  to  1. 
m  =  100 

DO  500  1  =  1,  m 
el  =  1 

lambda  =  0.01*el 

IF  (lambda  .gt.  0.9999)  lambda  =  0.9999 
x  =0.01 

c  Create  the  iteration  series  for  a  given  lambda 
N  =  7200 

WRITE  ( 1,  *) '  i  y ( i ) 

DO  200  i  =  1,  N 

y(i)  =  4. *lambda*x* ( 1-x) 
x  =  y( i) 

WRITE  (1,150)  i,  y ( i ) 

150  FORMAT  (lx,  14,  3x,  f8.5) 

200  CONTINUE 

c  With  iteration  series,  get  lengths  and  self-affine  fractal  dimension 
WRITE  (2,  *)'  Ln  (  e  )  Ln  (L(e))’ 
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c  Get  total  length  as  sum  of  differences  for  different  spacings  of  i 
DO  400  j  =  1,  48 
sum( j , 1)  =0. 

DO  300  i  =  3600  +  k(j),  7200,  k(j) 

sum( j , 1 )  =  sum ( j , 1 )  +  abs(  y(i)  -  y(  i-k(j)  )  ) 
d  IF  (j  .gt.  44) 

d  &  WRITE (6,250)  y ( i ) , y ( i-k( j ) ) , sum( j , 1 ) , i, j , k( j ) 

d250  FORMAT  (  lx,  3f8.4,3I5) 

300  CONTINUE 

d  IF  (j  .gt.  44) 

d  &  WRITE ( 6, * ) ' sum( ' , j , ' , 1 )  =  ’,sum(j,l) 

lglngth(j,l)  =  alogl0(  sum(j,l)  ) 
lgk(j)  =  alogl0(  float(k(j))  ) 
xy(l.j)  =  lgk ■ 4 ' 
xy(2,j)  =  lglngch(j,l) 

WRITE  (2,350)  lgk(j),  lglngth(j,l) 

350  FORMAT  (lx,  f8.5,  3x,  f8.5) 

400  CONTINUE 

c 

c  Do  linear  regression  to  get  Da  from  loglO(k)  vs.  loglO  (L(e)]  values 
c 

CALL  STLNRG ( n,  nyval,  xy,  ans) 
c 

c  -ans(2)  contains  Da 

c  ans (14)  contains  the  standard  deviation  of  Da 
c 

c  Put  lambda.  Da,  and  oDa  in  a  3  by  1  array 

c 

Dalmbda(l,l)  =  -ans (2) 

Dalmbda (2,1)  =  -ans (14) 

Dalmbda(3,l)  =  lambda 
500  CONTINUE 

WRITE  (3,550)  ( (Dalmbda(i, 1) ,  i  =  1,3),  1  =  l,m) 

550  FORMAT  (lx,  3(f9.5,2x)) 

STOP 

700  WRITE  (6,*) ’error  opening  feigenis.dat’ 

STOP 

800  WRITE  (6,*) ’error  opening  feigenll.dat’ 

STOP 

END 

c*********************************************************************** 

subroutine  stlnrg  (k,  nyval,  xy,  ans) 
c 

c  This  subroutine  computes  the  slope  and  other  statistics  for  a 
c  linear  regression  with  several  y  values  for  each  x  value  or  with 
c  one  independent  variable, 
c 


c 

argument 

use 

description 

c 

k 

input 

Number  of  different  x  values 

c 

c 

nyval 

input 

Array  of  number  of  y  values 

for 

each  x 

c 

value.  Dimensioned  k 

c 

xy 

input 

array  of  x  and  y  pairs  (xl, 

yi. 

x2 ,  y2 ,  etc 

c 

r* 

of  dimension  (2,k) 

c 

c 

ans 

output 

Array  of  results  computed,  dimensioned  16. 

c 

ans 

c 

1  Number 

c  *2  Slope 

c  3  Mean  of  y 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c* 


4  Mean  of  x 

5  y-intercept 

6  Sum  of  y**2 

7  Sum  of  squares  mean 

8  Sum  of  squares  slope 

9  Residual 

10  standard  deviation  of  x 

11  standard  deviation  of  y 

12  standard  deviation  error 

13  standard  deviation  y-bar 
*14  standard  deviation  slope 

15  standard  deviation  y-intercept 

16  F-Ratio  for  slope 


DIMENSION  nyval(l),  xy(l),  ans(16) 
sumx  -  0.0 
sumy  =  0.0 
sumx 2  =  0.0 
ans(6)  =  0.0 
sumxy  =  0.0 
n  =  0 
nx  =  1 
ny  =  2 

DO  20  j  =  l,k 
nyval(j)  =  1 
m  =  nyval ( j ) 
n  =  n  +  m 


DO  10  i  =  1,  m 

sumx  =  sumx  +  xy(nx) 

sumx2  =  sumx2  +  xy (nx) *xy (nx) 

sumy  *  sumy  +  xy(ny) 

ans(6)  =  ans(6)  +  xy(ny) *xy (ny) 

sumxy  =  sumxy  +  xy(nx) *xy (ny) 

10  ny  =  ny  +  1 

nx  =  nx  +  nyval(j)  +  1 
20  ny  =  ny  +  1 
en  =  n 
ans(l)  =  en 

si  =  ans(l)*sumx2  -  sumx*sumx 

s2  =  ans(l)*sumxy  -  sumx* sumy 

enl  =  en  -  1.0 

en2  =  enl  -  1.0 

ans(2)  =  s2/sl 

ans(3)  =  sumy/en 

ans(4)  =  sumx/en 

ans(5)  =  ans(3)  -  ans ( 2 ) *ans ( 4 ) 

ans(7)  =  ans(3)*sumy 

ans(8)  =  ans(2)*s2/en 

ans (9)  =  ans (6)  -  ans (7)  -  ans (8) 

s4  =  ans(9)/(en  -  2.0) 

endsl  =  en/sl 

ans(10)  =  sqrt( 1.0/ (enl*endsl) ) 

ans(ll)  =  sqrt((ans(6)  -  ans(7))/enl) 

ans (12)  =  sqrt(s4) 

al3  =  s4/en 

ans(13)  =  sqrt(al3) 

al4  =  s4*endsl 

ans(14)  =  sqrt(al4) 

ans(15)  =  sqrt(al3  +  al4*ans ( 4 ) *ans ( 4 )  ) 
ans(16)  =  ans(8)/s4 
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RETURN 

END 


noonnoo  ooo  no  non  o  ooo  ooo 


APPENDIX  B.  HENON  EQUATION  ANALYSIS  PROGRAM 


Q  **************************************************************** 

PROGRAM  CHAOS 


Written  by  Ray  Kamada  &  Korey  Jackson 

DOUBLE  PRECISION  X(5000),  Y(5000),  DA,  A,  B, 

+  DAERR,  L(100),  Bl,  S 

INTEGER  NMAX ,11, JMAX , IFLAG ,  Jl,  NMAX1 ,  COUNT.  EPSIL (100) 
OPEN ( 2 ,  FILE=' LENGTH. DAT ' ,  STATUS= 1  NEW ' ,  ERR=150) 

WRITE ( 2 , * ) '  A  B  DA  DAERR1 

OPEN ( 7 ,  FILE= 1 ENTROP . DAT 1 ,  STATUS= ' NEW’ ) 

WRITE ( 7 , * )  1  A  B  S’ 

Establish  NMAX 


NMAX=3600+1 
NMAX1=NMAX+100 
JMAX=NMAX/2 
PRINT* , 1 JMAX= 1 , JMAX 
IF  ( NMAX. GT. 5000)  GO  TO  100 
Initialize  all  variables 
I FLAG =0 

CALL  D I VI SR (NMAX,  JMAX,  EPSIL,  COUNT) 

Now  set  up  a  chaos  generating  routine... 

A=1 . 0570D0 
B1=0 . 3D0 
11=0 

0  IF ( II .GT. 600 ) GOTO  55 

J1=0 

A=A+0 . 0010D0 
B=B1 

11=11  +  1 

35  IF(J1.GT.30) GOTO  50 

J1=J1+1 
B=B+0 . 05D0 

print*, ’MAINl  A,  B, ' , A, B 
I FLAG =0 

CALL  HENON (NMAX1 , A, B, X, Y, IFLAG ) 
print* , ' hencn  complete ' , a, b, if lag 
Error  trap  for  x,  y  out  of  range  (diverging  solutions) 

I F ( I FLAG . EQ . 1 ) GOTO  35 

Comment  out  above  line  &  replace  w/below  line  when  B  is  fixed 
I F ( I FLAG . EQ . 1 ) GOTO  30 
if ( if lag. eq. 1 . ) goto  200 

Analyze  the  data 

CALL  ANAL Y S ( NMAX , JMAX , X , Y , A , B ,  DA,  DAERR,  COUNT,  EPSIL) 
print* analysis  complete',  a,  b,  DA 

44  WRITE (2, 45)  A,  B,  DA,  DAERR 

45  FORMAT (  F7 . 4 ,  3X,  F6.3,  5X,  F8.5,  5X,  F8.5) 

46  print* WRITE  to  file  complete  in  analysis' 

CALL  ENTROP ( NMAX, X,Y,S) 
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I 


WRITE;”, 47)  A,  B,  S 
47  FORMAT (  F7 . 4 , 3X,  F6 . 3 , 5X , F8 . 5 ) 

c  49  GOTO  35 
50  CONTINUE 
C  54  GOTO  30 
C 

55  CONTINUE 

C  CLOSE  (  2 ) 

C  CLOSE ( 7 ) 

GOTO  200 
C 

100  PRINT*, 'ERROR. .. .NMAX  Greater  than  10000 ... .Array  size  too  small' 
GOTO  200 

146  PRINT* ,’ ERROR. . .Can  not  open  henon.dat  file' 

GOTO  200 

148  PRINT* ,' ERROR. . .Can  not  open  analys.dat  file' 

GOTO  200 

150  PRINT*, 'ERROR. . .Can  not  open  length.dat  file' 

200  PRINT* ,' COMPLETED  MAIN2 ' 

END 

Q  ***************************************************************** 

SUBROUTINE  DIVISR ( NMAX , JMAX , EPSIL,  COUNT) 

INTEGER  CHECK,  CHECK1,  COUNT,  NMAX,  JMAX,  EPSIL(IOO) 

CHECK=0 
CHECK1=0 
COUNT=0 
MMAX=NMAX- 1 
1000  CHECK=CHECK+1 

IF ( CHECK. GT. JMAX)  GOTO  1050 

CHECK1=NMAX/CHECK 

CHECK1 =CHECK1 *CHECK 

IF ( CHECK1 . NE . NMAX)  GOTO  1000 

COUNT=COUNT+l 

EPSIL ( COUNT ) =CHECK 

PRINT*, 'EPSIL=' , EPSIL ( COUNT) ,  COUNT 
GOTO  1000 
1050  CONTINUE 

NMAX=NMAX+1 

RETURN 

END 

£  ***************************************************************** 
SUBROUTINE  HENON ( NMAX1 , A , B, X , Y , I FLAG) 

DOUBLE  PRECISION  X(NMAXl),  Y(NMAXl),  A,  B,  CHECK 
INTEGER  NMAX1,  N,  I FLAG,  NMAX 
C  OPEN ( 4, FILE=' HENON. DAT ' ,  STATUS= ’ NEW ’ ,  ERR=700) 

C 

C  Initialize  variables 
C 

PRINT*, 'PASSED  TO  HENON  OK’ 

DO  300  N=1 , NMAX1 
X ( N ) =0 . 0D0 
Y  ( N ) =0 . 0U0 
300  CONTINUE 
C 

C  MAIN  COMPUTATION  ROUTtne 
C 

C  First  initiate  X(l),  if  desired. 

X ( 1 ) =-0 . 65d0 
Y ( 1 ) =  0 . 38d0 
C 

c  WRITE (4  *) '  N  X  X ' 
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DO  400  N=1 ,  NMAX1 

X ( N+l ) =1 . -A*X (N ) *X (N) +Y ( N ) 

CHECK=X ( N+l ) 

IF ( CHECK. GE. 100. OR. CHECK. LE. -100. )  THEN 
I FLAG= 1 

PRINT*, 'PROGRAM  OUT  OF  RANGE  IN  X’ 

END  IF 

Y (N+l ) =B*X ( N) 

CHECK=Y(N+1) 

IF(CHECK.GE. 1000. OR. CHECK. LE. -1000.0)  THEN 
IFLAG=1 

PRINT* ,  *  PROGRAM  OUT  OF  RANGE  IN  Y ’ 

END  IF 

IF ( IFLAG . EQ. 1 ) GOTO  760 
IF(N.LT. 101)GOT0400 
C  WRITE (4,350)  N,  X(N),  Y (N) 

350  FORMAT (I10,5X,E15.5,5X,E15.5) 

400  CONTINUE 
C 

C  NOW  SHIFT  THE  VALUES  DOWN  IN  X(Y),  Y(I) 

C 

NMAX=NMAX1-100 
DO  600  N=1 ,  NMAX 
X ( N ) =X ( N+100 ) 

Y ( N ) =Y (N+100 ) 

600  CONTINUE 
C 

C  Return  to  the  main  program 
C 

PRINT*, 'COMPLETED  HENON  OK  FOR’,A,B 
GO  TO  770 
C 

700  PRINT*, ’Error  in  opening  HENON . DAT  file" 

C  740  CLOSE (2) 

C  750  CONTINUE 

GOTO  770 

760  PRINT*, 'X  value  out  of  range... A, B  TOSSED’,  A,  B 
IFLAG=1 

770  CLOSE ( 1 ) 

RETURN 

END 

SUBROUTINE  ANALYS ( NMAX , JMAX ,  X,Y,A,B,  DA,  DAERR,  COUNT, 

+  EPSIL) 

DOUBLE  PRECISION  SCRAP1,  SCRAP2 ,  L(100),  X(NMAX),  SUM(100), 
+  Y(NMAX),  EPS ( 100 ) ,  DA,  A,  B,  DAERR,  SCRAP3 ,  XY(200), 

+  ANS ( 16 ) 

INTEGER  NMAX,  I,  J,  CHECK,  CHECK1 , JMAX , COUNT,  EPSIL(IOO), 

+  COUNT2 ,  nyval ( 100 ) ,  COUNT3 
C 

C  VARIABLE  LIST: 

C  NMAX=MAXIMUM  NUMBER  OF  DATA  POINTS 

C  COUNT=NUMBER  OF  EVEN  DATA  POINTS  THAT  WILL 

C  DIVIDE  EVENLY  INTO  NMAX 

C  L(J)=LENGTH  FOR  A  GIVEN  EPSILON 

C  EPSIL (J ) =WIDTH,  OR  AN  INTEGER  DIVISOR  OF  NMAX 

C  EPS( J)=EPSIL(J)  IN  INTEGER  FORM 

C  DAERR=STANDARD  DEVIATION  OF  DA  (LOG  FORM) 

C 

OPEN  (3,  FILE= 'ANALYS. DAT ’ ,  STATUS= ' NEW ’ , 

+  ERR=90 10 ) 


128 


c 

C  INITITALIZE  THE  VARIABLES 

C 

DO  4000  J=l,  100 
L( J)=0.0D0 
SUM( J) =0.0D0 
4000  CONTINUE 
C 

C  Now  compute  L(epsil) 

C 

DO  7500  J=l,  COUNT 
CHECK=NMAX -E P S I L { J ) 

CHECK1=EPSIL(J) 

DO  7000  I=CHECK1,  CHECK,  CHECK1 
SCRAP 1=X ( 1  +  1 ) -X  < I-CHECK1+1 ) 

SCRAP2=Y ( 1+1 ) -Y ( I-CHECK1+1 ) 

SCRAP1=SCRAP1*SCRAP1 
SCRAP2=SCRAP2*SCRAP2 
SCRAP2=SCRAP2+SCRAP1 
SCRAP2=DSQRT ( SCRAP2 ) 

SUM ( J ) =SUM ( J ) +SCRAP2 
7000  CONTINUE 

L( J. =SUM( J) 

7500  CONTINUE 
C 

C  QUIK  DATA  PRINT 

DO  7550  J=l,  COUNT 
L(J)=L( J)+1.0D0 
SCRAP1=L(J) 

L(J)=  DLOG10 ( SCRAP1 ) 

SCRAP3=DBLE ( EPSIL ( J ) ) 

EPS ( J) =DLOG10 ( SCRAP3 ) 

WRITE (3,7545)  EPS(J),  L(J) 

7545  FORMAT ( 5X , F8 . 5 , 5X , F8 . 5 ) 

7550  CONTINUE 
C 

C  Now  do  a  linear  regression  to  find  DA 

C  First  put  into  a  format  to  use  the  canned  linear  regression 
C  subroutine. 

COUNT3=2  *COUNT 
DO  7700  J=2 ,  COUNT 3 ,  2 
I=J/2 

XY ( J-l ) =EPS ( I ) 

XY ( J) =L( I ) 

7700  CONTINUE 
C 

C  COUNT2  IS  THE  COUNT  WITH  THE  LAST  3  VALUES  TOSSED  FOR  Da 
c  computation 

COUNT2=COUNT-3 

CALL  LINREG  ( COUNT2 ,  nyval,  XY,  ANS ) 
print*, 'made  it  out  of  linreg' 

DA=ANS ( 2 ) 

DA=-DA 

DAERR=ANS( 14) 

GOTO  9100 

9000  PRINT*, 'ERROR  OPENING  HENON.DAT  FILE' 

9010  PRINT*, 'ERROR  OPENING  ANALYS.DAT  FILE' 

9020  PRINT*, 'ERROR  OPENING  LENGTH . DAT  FILE’ 

9100  print*, 'made  it  to  9100' 

RETURN 

END 
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subroutine  LINREG  (k,  nyval,  xy,  ans) 


c 

c 

This  subroutine  computes  the  slope  and  other  statistics  for  a 

c 

linear  regression 

with  several  y  values  for  each  x 

value  or  with 

c 

r* 

one  independent  variable. 

c 

argument 

use 

description 

c 

k 

input 

Number  of  different  x  values 

c 

c 

nyval 

input 

Array  of  number  of  y  values  for 

each  x 

c 

value.  Dimensioned  k 

c 

xy 

input 

array  of  x  and  y  pairs  (xl,  yl. 

x2 ,  y2 ,  etc . ) 

c 

of  dimension  (2,k) 

c 

c 

c 

ans 

output 

Array  of  results  computed,  dimensioned  16. 

c 

ans 

c 

1  Number 

c 

*2  Slope 

c 

3  Mean  of  y 

c 

4  Mean  of  x 

c 

5  y-intercept 

c 

6  Sum  of  y**2 

c 

7  Sum  of  squares  mean 

c 

8  Sum  of  squares  slope 

c 

9  Residual 

c 

10  standard  deviation  of  x 

c 

11  standard  deviation  of  y 

c 

12  standard  deviation  error 

c 

13  standard  deviation  y-bar 

c 

*14  standard  deviation  slope 

c  15  standard  deviation  y-intercept 

c  16  F-Ratio  for  slope 

c  **************************************************************** 

DOUBLE  PRECISION  nyval (1),  xy(l),  ans (16),  si,  s2,  sumx,  sumx2 
+  sumy,  sumy2,  sumxy,  enl,  en2,  s4,  endsl,  al3,  al4 
sumx  =  0.0D0 
sumy  =  0.0D0 
sumx2  =  0.0D0 
ans (6)  =  0.0D0 
sumxy  =  0.0D0 
n  =  0 
nx  =  1 
ny  =  2 

DO  20  j  =  1, k 

nyval ( j )  =  1 

m  =  nyval ( j ) 

n  =  n  +  m 

DO  10  i  =  1,  m 

sumx  =  sumx  +  xy(nx) 

sumx2  =  sumx2  +  xy ( nx) *xy ( nx ) 

sumy  =  sumy  +  xy ( ny ) 

ans(6)  =  ans(6)  +  xy ( ny ) *xy ( ny ) 

sumxy  =  sumxy  +  xy(nx) *xy(ny) 

10  ny  =  ny  +  1 

nx  =  nx  +  nyval(j)  +  1 
20  ny  =  ny  +  1 
en  =  n 
ans(l)  =  en 

si  =  ans  ( 1 )  *surnx2  -  sumx*sumx 
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92  =  ang ( 1 ) *aumxy  -  aumx*aumy 

enl  =  en  -  1.0 

en2  =  enl  -  1.0 

ans(2)  =  s2/sl 

ana(3)  =  sumy/en 

ana (4)  =  sumx/en 

ans(5)  =  ana (3)  -  ana (2 ) *ana ( 4) 

ana(7)  =  ana(3)*aumy 

ana(8)  =  ana (2 ) *a2/en 

ana (9)  =  ana (6)  -  ana(7)  -  ana (8) 

a4  =  ana(9)/(en  -  2.0) 

endal  =  en/al 

ana(lO)  =  sqrt ( 1 . 0/ (enl*endal ) ) 

ana (11)  =  aqrt((ana(6)  -  ana(7))/enl) 

ana (12)  =  aqrt(a4) 

al3  =  a4/en 

ana (13)  =  sqrt(al3) 

al4  =  a4*endal 

ana (14)  =  aqrt(al4) 

ana(15)  =  aqrt(al3  +  al4*ana (4 ) *ana ( 4)  ) 

ana(16)  =  ana(8)/a4 

RETURN 

END 

Q  *****★**********★★★****★★★★★★★★***★★*★**★★★★★★*********■**★★**★ 

C  23  Apr  92  reviaion  (goea  with  0<a<1.25  plot) 

C 

SUBROUTINE  ENTROP ( NMAX ,  X ,  Y ,  S  ) 

DOUBLE  PRECISION  X(NMAX),  Y (NMAX) ,  S,  XMAX,  XMIN,  YMAX,  YMIN, 
+  PROB ( 200 , 200 ) ,  INVNMA,  P,  X RANGE ,  Y RANGE 
INTEGER  IX,  IY,  I,  J,  NMAX 
C  Firat  find  maximum  and  minimum  valuea 
XMAX=0 . 0D0 
YMAX=0 . 0D0 
XMIN=0 . 0D0 
YMIN=0 . 0D0 

print* finding  max  and  minimum  valuea' 

DO  2000  1=1,  NMAX 

IF (X ( I ) .GT. XMAX)  XMAX=X ( I ) 

IF (X ( I ) . LT. XMIN)  XMIN=X( I ) 

IF(Y(I) .GT.YMAX)  YMAX=Y(I) 

IF ( Y ( I ) . LT. YMIN)  YMIN=Y(I) 

2000  CONTINUE 

C  Now  initialize  the  probability  array 
print* ,' initializing  P  array' 

DO  2105  1=1,  200 
DO  2100  J=l,  200 
PROB ( I , J) =0 . 0D0 
2100  CONTINUE 

2105  CONTINUE 

C  Count  the  number  of  pointa  in  each  unit  area 
C  Firat  aet  up  the  acheme 
XRANGE=XMAX-XMIN 
XRANGE=2 . 00D2/XRANGE 
YRANGE= YMAX- YMIN 
YRANGE=2 . 00D2/YRANGE 
DO  2110  1=1,  NMAX 

IX=INT ( ( X ( I ) -XMIN ) *XRANGE ) + 1 

IY=INT ( (Y(I)-YMIN) *YRANGE ) +1 

IF ( IX . GT . 200 )  IX=200 

IF ( IY . GT. 200 )  IY=200 

PROB ( Iv  IY)=PROB(IX,IY)+1.0D0 


131 


2110  CONTINUE 
C  Now  normalize 

INVNMA=1 . 0D00/DBLE ( NMAX ) 

DO  2120  IX=1,  200 
DO  2115  IY=1 ,  200 

PROB ( IX,  I Y ) =PROB (IX, IY ) * INVNMA 
2115  CONTINUE 

2120  CONTINUE 

C  Finally,  compute  the  entropy  S 
S=0 . 0D0 

print* computing  S' 

DO  2130  IX=1,  200 
DO  2125  IY=1,  200 

IF ( P . ne . OdO ) print* , 'P=' ,P 
P=PROB ( IX, IY ) 

IF(P.EQ.ODO)  GOTO  2125 
S=s-P*DLOG(P) 

2125  CONTINUE 

2130  CONTINUE 

C  Now  return  to  the  main  program  to  write  S  for  this  value  of  A 
RETURN 
END 
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APPENDIX  C.  PROGRAM  MCNIDER 


PROGRAM  MCNIDE 

C 

c  This  program 

c  1)  computes  Monte  Carlo  particle  velocity  and  position  as  a  function 
c  of  time  step,  t+A*t,  using  the  McNider  algorithm,  Kamada  vertical 

c  velocity  skewness  variations,  and  Kamada  mesoscale  windflow  and 

c  turbulence  simulation  model, 

c 

c  2)  computes  the  chaos  metrics.  Da  (self-affine  fractal  dimension) 

c  for  velocity,  S  (information  entropy),  and  lambda  (the  Lyapunov 

c  exponent),  as  well  as  Ri,  the  atmospheric  Richardson  number,  BVF 

c  ( Brunt -Vaisala  frequency)  1  (buoyancy  length  scale),  ou,v,w,  tne 

c  velocity  variances,  TKE ,  the  turbulence  kinetic  energy,  Db 

c  (self-similar  fractal  dimension  for  position), 

c 


c 

DIMENSION  nyval ( 45 ) ,  k<45),  ipts(lOO),  bi(2,1800) 

DOUBLE  PRECISION  Linv,  lglngth(45,  100),  lgk(45), 

&  vert ( 0 : 3700 ) ,  r(0:3700),  rl(0:3700),  r2(0:3700), 

&  el,  p(1000),  entropy(lOl) ,  up,  dn,  lyap,  pu,  pv,  pw, 

&  sum (45,  100),  DavsL(8, 100) ,  pk,  ans(16),  xy(90) 

EQUIVALENCE  (vert,  r  ) 

DATA  k/1,  2,  3,  4,  5,  6,  8,  9,  10,  12,  15,  16,  18,  20,  24,  25, 

&  30,  36,  40,  45,  48,  50,  60,  72,  75,  80,  90,  100, 

&  120,  144,  150,  180,  200,  225,  240,  300,  360,  400, 

&  450,  600,  720,  900,  1200,  1800,  3600/ 

OPEN  (1,  FILE= ' McNideis . dat ' ,  RECL=255,  STATUS= ' new • ,  ERR=700) 

OPEN  (2,  FILE=’ McNidell.dat' ,  RECL=255,  STATUS= ’ new ' ,  ERR=800) 

OPEN  (3,  FILE= ' DavsL . dat ' ,  RECL=255,  STATUS= ’ new ' ,  ERR=900) 

OPENa  (4,  FILE= ' SvsL . dat ' ,  RECL=255,  STATUS= ’ new ’ ,  ERR=1000) 

OPEN  (7,  FILE= ' lyap . dat ' ,  RECL=255,  STATUS= ' new 1 ,  ERR=1100) 

OPEN  (8,  FILE=' bifurct.dat' ,  RECL=2o5,  STATUS= ' new ' ,  ERR=1200) 

OPEN  (9,  FILE= ' initial . dat ' ,  RECL=255,  STATUS^ ' old ' ,  ERR=1300) 
c 

c  INITIALIZE 

q  *  *  *  *  -k  ★  *  *  ★  * 

READ  (9,30) 

30  FORMAT  (///) 

READ  (9,*)  iis,  ill,  iDavsL,  isvsL,  ilyap,  ibifurct,  noskew, 

St  iMcNid,  iKamada,  dt,  ivert,  iphase,  iveloc,  windspd2,  zO, 

St  Start,  m2,  iw,  ix,  upfactor,  dnfactor,  dwfactor,  zinitial, 

St  icount,  iposi 

WRITE ( 6 , * ) ' iis ,  ill,  iDavsL,  isvsL,  ilyap,  ibifurct,  noskew,' 

WRITE ( 6 ,*)' iMcNid,  iKamada, dt, ivert, iphase, iveloc,  windspd2 , zO , ' 
WRITE ( 6, *)' Start ,  m2,  iw,  ix, upfactor,  dnfactor, dwfactor, zinitial ' 
WRITE ( 6, *)' icount,  iposi' 

WRITE ( 6 , * )  iis,  ill,  iDavsL,  isvsL,  ilyap,  ibifurct,  noskew, 

St  iMcNid,  iKamada,  dt,  ivert,  iphase,  iveloc,  windspd2,  zO, 

St  Start,  m2,  iw,  ix,  upfactor,  dnfactor,  dwfactor,  zinitial, 

St  icount,  iposi 

c  Comments  on  input  parameters: 

c  iis  =  1  switch  gives  position  and/or  velocity  file 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


iil  =  1  switch  gives  Ln  e  vs  Ln(L(e)  file 

iDavsL  =  1  switch  gives  Da  vs  1/L  file 

isvsL  =  1  switch  gives  entropy,  S,  vs  1/L  file 
ilyap  =  1  switch  gives  Lyapunov  exponent  vs  1/L  file 

ibifurct  =  1  switch  gives  bifurcation  disagram  file 

noskew  =  1  switch  turns  off  vertical  velocity  skewness 
computation 

iMcNid  =  1  switch  turns  on  McNider  vertical  velocity  skewness 
scheme 

iKamada  =  1  switch  turns  on  Kamada  vertical  velocity  skewness 
scheme 

dt  =  time  step  size  in  seconds 

ivert  =  1  switch  gives  vertical  height  file 

iphase  =  1  switch  gives  phase  space  vector  file 

iveloc  =  1  switch  gives  velocity  vector  file 

windspd2  =  input  proscribed  wind  speed  at  2  meters 

zi  =  input  proscribed  fully  turbulent  boundary  layer 

height 

zO  =  input  proscribed  surface  roughness  height  (“1/7)  the 

average  height  of  surface  roughness  elements  if 
element 


density  is  between  10  and  50%  (like  most  vegetation) 


Do  some  idiot  proofing 


IF  (  ivert  .eq.  1)  THEN 
iphase  =  0 
iveloc  =  0 
iposi  =  0 
ENDIF 

IF  (  iphase  .eq.  1)  THEN 
ivert  =  0 
iveloc  =  0 
iposi  =  0 
ENDIF 

IF  (  iposi  .eq.  1)  THEN 
ivert  =  0 
iphase  =  0 
iveloc  =  0 
ENDIF 

IF  (  noskew  .eq.  1)  THEN 
iMcNid  =  0 
iKamada=  0 
ENDIF 

IF  (  iMcNid  .eq.  1)  THEN 
noskew  =  0 
iKamada=  0 
ENDIF 

IF  (  iKamada. eq.  1)  THEN 
ncskew  =  0 
iMcNid  =  0 
ENDIF 

Initialize  variables 
ml  =  1 
n  =  3600 

en  =  n 

WRITE  (4, * )  ’  S  1/L' 

WRITE  (7,*)'  Lyap  1/L' 

WRITE  (8,*)  •  1/L  bifurcation  pt ' 

DO  500  1  =  ml,  m2 

Let  Linv,  the  inverse  Obukhov  length,  1/L  =  -kgw'O’ / (0*ustar**3 ) ,  be 
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ir>  m  m  sn 


c  varied  in  increments  along  the  interval  from  -1  to  1. 
em2  =  m2 

Linv  =  start  +  0. 3d0* ( 1-1 ) /em2 

zi  =  2000  -  6000* (Linv  +  .2) 

r(0)  =  0. 

rl(0)  =  0. 

r2 (0)  =  0. 

vert(0)  =  zinitial 

x  =  0.001 

y  =  0.002 

z  =  zinitial 

udp  =  0.001 

vdp  =  0.002 

wdp  =  0.003 

Q  ******** 

c  GET  iteration  series  for  a  given  1/L 
0  ******** 

IF  (  iis  .ne.  1  )  GOTO  150 
pu  =  0.0 
pv  »  0.001 

c  pw  =  -.01*alogl0( 1) 

pw  =  -0.001 
time  =  0.0 
rmax  =  0.0 
r2max  =  0.0 
rmin  =  0.0 
r2min  =  0.0 
sigU2sum  =  0. 
sigV2sum  =  0. 
sigW2sum  =  0. 
zsum  =  0. 

BVFsum  *  0. 

BLSsum  =  0. 
sumdw  =  0. 
tkebar  =  0. 
wbuoy  =  0. 
sumw  =  0 . 
sumwbuoy  =  0 . 

d  WRITE ( 6, * ) ' rmax 2  =',rmax2 

rmax2  =  0 . 
rmin2  =  0. 
iflag  =  0 

DO  100  i  =  1,  n  +  100 

CALL  McNid  (x,  y,  z,  udp,  vdp,  wdp,  dt,  noskew,  iMcNid, 
iKamada,  Linv,  windspd2,  deltaz,  zi,  zO,  wbuoy,  sumw,  Ri, 
sumwbuoy,  BVF,  BLS,  sigmau2,  sigmav2,  sigmaw2,  tke,  icount, 
pu,  pv,  pw,  i,  time,  wk,  ustar,  iw,  ix,  upfactor,  iflag, 
sigmau,  sigmav,  sigmaw,  Tlw,  sumdw,  dnfactor,  dwfactor  ) 
BVFsum  =  BVFsum  +  BVF 
BLSsum  =  BLSsum  +  BLS 
sigU2sum  =  sigU2sum  +  sigmau**2 

sigV2sum  =  sigV2sum  +  sigmav**2 

sigW2sum  =  sigW2sum  +  sigmaw**2 

zsum  =  zsum  +  z 

bltime  =  0. 3*zi/ ( sigmaw*3700 . ) 
eddytime  =  0.01*Tlw 
c  IF  (Linv  .ge.  0.000)  dt  =  MIN  (eddytime,  bltime) 

c  IF  (Linv  .ge.  0.000)  dt  =  MIN  (dt,  0.25) 

c  IF  (Linv  . gt .  0.000)  dt  =  eddytime 

c  IF  (Linv  .It.  -0.00)  dt  =  MIN  ( 1 . 2*zi/ ( 3700 . 0*wk) ,  eddytime) 

dt  =  0.166666667 
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time  =  time  +  dt 

cd  WRITE  (6,*)  'time', time 

c  get  position  vector,  rl 

rl(i)  =  sqrt(  x*x  +  y*y  +  ( z-zinitial ) **2  ) 
c  get  velocity  vector,  r2 

r2(i)  =  sqrt(  udp*dup  +  vdp*vdp  +  wdp*wdp  ) 
c  get  vertical  position  series 

IF  (  ivert  .eq.  1)  r(i)  =  z  -  zinitial 
c  get  position  series 

IF  (  iposi  .eq.  1)  r(i)  =  rl(i) 
c  get  phase  space  position 

IF  (iphase.eq.  1)  r(i)  =  sqrt(  rl(i)*rl(i)  +  r2(i)*r2(i)  ) 
c  get  velocity  series 

IF  (iveloc  .eq.  1)  r(i)  =  r2(i) 
c  get  min/max  range  of  the  values 

IF  (  i  .gt.  100  .and.  r(i)  .gt.  rmax  )  rmax  =  r(i) 

IF  (  i  .gt.  100  .and.  r(i)  .It.  rmin  )  rmin  =  r(i) 

IF  (  i  .gt.  100  .and.  r(i)  .gt.  rmax2)  rmax2=  r(i) 

IF  (  i  .gt.  100  .and.  r(i)  .It.  rmin2)  rmin2=  r(i) 

c  WRITE  (1,50)  i,  rl(i),  r2(i),  r(i) 

c50  FORMAT  (lx,  14,  3(3x,  fl2.5)  ) 

100  CONTINUE 

150  CONTINUE 

WRITE  (6, *)' iteration  series  created  for  1/L  =  ',  Linv 
IF  (  ibifurct  .ne.  1  )  GOTO  240 

Q  *************** 

c  GET  BIFURCATION  diagram,  using  the  last  1600  points.  For  each  1/L 
value 

c  put  all  the  unique  points  in  "bi" 

Q  *************** 

C 

m  =  1 

bi(l,m)  =  Linv 
bi ( 2 , m)  =  r ( 1999 ) 
slop  =  0.01* (rmax  -  rmin) 
c  IF  (iveloc  .eq.  1)  slop  =  0.005*11. 

c  IF  (ivert  .eq.  1)  slop  =  0.005*zi 

c  IF  ( iphase  .eq.  1)  slop  =  0.002*3000. 

c  Check  every  r(i) 

DO  220  i  =  2000,  n+99 

c  Check  each  r(i+l)  against  all  unique  points  in  bi 

DO  180  j  =  1,  m 

up  =  bi(2,j)  +  slop 
dn  =  bi ( 2 , j )  -  slop 

IF  (  r(i+l)  .gt.  dn  .and.  r(i+l)  .It.  up  )  GOTO  220 
c  Since  r(i+l)  lies  beyond  margins  of  all  "bi",  we  assume  it  a  new 

c  unique  point,  so  add  it  to  the  list  of  values  in  "bi". 

180  CONTINUE 

bi(2,m+l)  =  r(i+l) 
bi(l,m+l)  =  Linv 
m  =  m+1 

220  CONTINUE 

c  IF  (  iveloc  .ne.  1  )  bi(2,m)  =  r(1999) 

WRITE  (8,230)  ((  bi ( i, j ) ,  i  =  1,2),  j  =  l,m) 

230  FORMAT  (  lx,  f9.5,  5x,  f9.5  ) 

240  CONTINUE 

WRITE  (6, *) 'bifurcation  diagram  points  found  for  1/L  =  ',  Linv 
IF  (  ilyap  .ne.  1  )  GOTO  260 

Q  ************ 

c  GET  LYAPUNOV  exponent  for  this  series, 

Q  ************ 
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I 

I 


,  where  f’(x  )  =/j(1  -2x  ) 

i+1  i+1 


c  n 

c  lyap  =  lim  (1/n)  sum  log  jf'(x  ) 

c  nxo  i=l  e  i+l 

c 

c  where  we  start  from  y(100)  to  avoid  dependence  on  initial  condition, 
c  xO. 

c 

lyap  =  0.0 
DO  250  i  =  1,  n 

c  lyap  =  lyap  +  dlog(  abs(mu*(1.0  -  2.*y(i+100)  )  )) 

lyap  =  0.0000000000 
250  CONTINUE 

c  Normalize  the  lyapunov  exponent  by  N  -  100,  the  number  of  points 
c  scanned. 

lyap  =  lyap/en 
WRITE  (7,255)  lyap,  Linv 
255  FORMAT  (lx,  2(f9.5,2x)  ) 

260  CONTINUE 

IF  (  isvsL  .ne.  1  )  GOTO  290 
c 

c  ***********  jj 

c  GET  ENTROPY  value,  S  =  -  sum  p  log  p  ,  for  this  iteration  series  &  1/L 

Q  ***********  j.  0  4 

c 

c  Zero  out  the  number  of  points  in  each  bin 
DO  265  ii  =  1,  100 
ipts(ii)  =  0 
265  CONTINUE 

c  Subdivide  the  interval  into  100  sub-intervals, 

c  ii,  and  check  all  pts(i),  from  101  -  3600. 

entropy (1)  =  0.0 
DO  270  i  =  101,  n  +  100 

ii  =  99* (  r(i)  -  rmin2 ) / (rmax2  -  rmin2)  +  1.01 
c  Make  sure  we  don't  create  extra  sub-intervals 

IF  (ii  .gt.  100)  ii  =  100 

c  Increment  the  ipts  counter  for  sub-interval  ii  for  every  r(i)  found 

c  in  ii. 

ipts(ii)  =  ipts(ii)  +  1 
270  CONTINUE 

DO  285  ii  =  1,  100 

c  If  no  points  are  in  sub-interval,  ii,  leave  the  entropy  unchanged. 

IF  (  ipts(ii)  .eq.  0)  THEN 
si  =  0.0 
GOTO  280 
ENDIF 

c  Otherwise  compute  probability  of  a  point  being  in  ii  as  the  #  of 

c  points  in  ii  divided  by  N,  the  total  number  of  points  scanned, 

pts  =  ipts(ii) 
prob  =  pts/en 
si  =  prob*alog(prob) 

c  Since  the  probability  is  less  than  or  equal  to  1,  the  log  is 

c  negative,  so  si  is  less  than  or  equal  to  zero,  so  the  entropy  will 

c  grow  more  positive,  if  the  r(i)  are  scattered  over  more  than  one 

c  sub-interval. 

280  entropy(l)  =  entropy(l)  -  si 

285  CONTINUE 

WRITE  (4,255)  entropy(l),  Linv 

IF  (  ill.  eq.  1)  WRITE  (2,  *)'  Ln(e)  Ln  (L(e)]’ 

WRITE  ( 6 ,*)' entropy  value  found  for  1/L  =  ',  Linv 
290  CONTINUE 

IF  (  iDavsL  .ne.  1  )  GOTO  500 
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c  GET  Da.  First  get  total  length  as  sum  of  y(i)  changes  for  varying 
c  *******  spacings  of  i. 
nn  =  45 

DO  400  j  =  1,  nn 

sum(j,l)  =  1.0000 

DO  300  i  =  k(j),  3600,  k(j) 

sum(j,l)  =  sum(j,l)  +  abs(  r(i+100)  -  r(i-k(j)+100  )) 

300  CONTINUE 

c  Get  log  of  the  total  length 

lglngth{j,l)  =  dlogl0(  aum(j,l)  ) 
c  Get  log  of  the  spacing 
pk  =  k(j) 

lgk(j)  =  dlogl0(  pk  ) 

c  Put  these  log  values  into  an  array  and  send  to  the  linear  regression 

c  routine. 

xy(2*j-l)  =  lgk(j) 
xy(2* j )  =  lglngth( j , 1 ) 

IF  (ill  .eq.  1)  WRITE  (2,350)  lgk(j),  lglngth(j,l) 

350  FORMAT  (lx,  f8.5,  3x,  f8.5) 

400  CONTINUE 

c  WRITE  (6,*)  'program  passed  statement  400' 

c 

c  DO  LINEAR  REGRESSION  to  get  Da  from  loglO(k)  vs.  loglO  [L(e)]  values 
c 

c  Use  only  spacings  up  to  1200  points  for  the  linear  regression.  So 
c  skip  spacings  1800  and  3600. 

nn  =  nn  -  2 

CALL  STLNRG ( nn ,  nyval,  xy,  ans) 


-ans(2)  contains  Da 
c  ans (14)  contains  the  standard  deviation  of  Da 
c 

c  Put  1/L,  Da,  oDa,  BVF,  BLS,  aw3,  tke  in  a  7  by  1  array 
c 


500 


550 


ptsinv  =  l./(n+100.) 

BVFbar  =  BVFsum*ptsinv 
BLSbar  =  BLSsum*ptsinv 
sigw2bar  =  sigw2 sum* ptsinv 

tkebar  =  ,5*(sigu2sum  +  sigv2sum  +  sigw2sum) *ptsinv 
zbar  =  zsun^ptsin” 

DavsL(l,l)  =  -ans(2) 

DavsL(2,l)  *  -ans (14) 

DavsL(3,l)  =  Linv 
DavsL(4,l)  =  BVFbar 
DavsL(5,l)  =  BLSbar 
DavsL(6,l)  =  sigw2bar 
DavsL(7,l)  =  tkebar 
DavsL(8,l)  =  zbar 

WRITE  (6,*) 'Da  value  found  for  1/L  =  ’,  Linv 
CONTINUE 


WRITE 

& 

(3,*)  ' 

1 

9 

Da 

ow2 

oDa 

tke 

Linv 
zavg ' 

BVF 

BLS’ 

WRITE 

& 

(6,*)  ' 

f 

9 

Da 

ow2 

a  Da 
tke 

Linv 
zavg ' 

BVF 

BLS’ 

WRITE  (3,550)  ( (DavsL( i, 1) ,  i  =  1,8),  1  =  ml, m2) 
WRITE  (6,550)  ( ( DavsL ( i , 1 ) ,  i  =  1,8),  1  =  ml, m2) 
FORMAT  (lx,  5 ( f 10. 4, lx) ,  2f7.2,  f7.0) 

STOP 

WRITE  (6,*) 'error  opening  McNideis.dat' 

STOP 


700 


800 

WRITE 

STOP 

(6,* 

900 

WRITE 

STOP 

(6,* 

1000 

WRITE 

STOP 

(6,* 

1100 

WRITE 

STOP 

<6,* 

1200 

WRITE 

STOP 

(6,* 

1300 

WRITE 

STOP 

END 

(6,* 

)'error  opening  McNidell.dat' 
) 'error  opening  DavsL.dat' 

) 'error  opening  SvsL.dat' 

) 'error  opening  Lyap.dat' 

)’ error  opening  bifurct.dat’ 

) 'error  opening  initial.dat’ 


c 

c 

c 

c 

c 


SUBROUTINE  STLNRG  (k,  nyval,  xy,  ana) 

This  subroutine  computes  the  slope  and  other  statistics  for  a 
linear  regression  with  several  y  values  for  each  x  value  or  with 
one  independent  variable. 


c 

argument 

use 

description 

c 

k 

input 

Number  of  different  x  values 

c 

c 

nyval 

input 

Array  of  number  of  y  values  for  each  x 

c 

value.  Dimensioned  k 

c 

xy 

input 

array  of  x  and  y  pairs  (xl,  yl,  x2,  y2 

c 

of  dimension  (2,k) 

c 

c 

c 

ans 

output 

Array  of  results  computed,  dimensioned 

c 

ans 

c 

1  Number 

c 

*2  Slope 

c 

3  Mean  of  y 

c 

4  Mean  of  x 

c 

5  y-intercept 

c 

6  Sum  of  y**2 

c 

7  Sum  of  squares  mean 

c 

8  Sum  of  squares  slope 

c 

9  Residual 

c 

10  standard  deviation  of  x 

c 

11  standard  deviation  of  y 

c 

12  standard  deviation  error 

c 

13  standard  deviation  y-bar 

c 

*14  standard  deviation  slope 

c 

15  standard  deviation  y-intercept 

c 

16  F-Ratio  for  slope 

V— 

DIMENSION  nyval(l) 

DOUBLE  PRECISION  xy(l),  ans(16) 
c  WRITE  (6,*)'  program  entered  linear  regression  routine' 

sumx  =  0.0 
sumy  =  0.0 
sumx2  =  0.0 
ans(6)  =  0.0 
sumxy  =  0.0 
n  =  0 
nx  =  1 
ny  =  2 

DO  20  j  =  1 , k 
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nyval(j)  =  1 

m  =  nyval(j) 

n  =  n  +  m 

DO  10  i  =  1,  m 

sumx  =  sumx  +  xy ( nx ) 

sumx2  =  sumx2  +  xy(nx) *xy (nx) 

aumy  =  sumy  +  xy(ny) 

ans(6)  =  ans(6)  +  xy (ny) *xy (ny) 

sumxy  =  sumxy  +  xy(nx) *xy(ny) 

10  ny  =  ny  +  1 

nx  =  nx  +  nyval ( j )  +1 
20  ny  =  ny  +  1 
en  =  n 
ans(l)  =  en 

si  =  ans(l)*sumx2  -  sumx*sumx 

s2  =  ans(l)*sumxy  -  sumx*sumy 

enl  =  en  -  1.0 

•en2  =  enl  -  1.0 

ans(2)  =  s2/sl 

ana (3)  =  aumy/en 

ana (4)  =  sumx/en 

ana(5)  =  ana (3)  -  ans(2)*ans(4) 

ans(7)  =  ana ( 3 ) *  sumy 

ans(8)  =  ans(2)*s2/en 

ana (9)  =  ana (6)  -  ans(7)  -  ans(8) 

s4  =  ans(9)/(en  -  2.0) 

endal  =  en/al 

an3(10)  =  sqrt ( 1 .0/ (enl*endsl) ) 

ana (11)  =  sqrt ((ana (6)  -  ana(7))/enl) 

ana (12)  =  sqrt (34) 

al3  =  s4/en 

ana (13)  =  sqrt(al3) 

al4  =  s4*endsl 

ana (14)  =  aqrt(al4) 

ans(15)  =  sqrt(al3  +  al4*ans ( 4 ) *ans ( 4 )  ) 

ans(16)  =  ana(8)/a4 

RETURN 

END 

SUBROUTINE  McNid  (x,  y,  z,  udp,  vdp,  wdp,  dt,  noskew,  iMcNid, 

&  IKamada,  Linv,  windapd2,  deltaz,  zi,  zO,  wbuoy,  sumw,  Ri, 

&  sumwbuoy,  BVF,  BLS,  sigmau2,  aigmav2,  aigmaw2,  tke,  icount, 

&  pu,  pv,  pw,  i,  time,  wk,  ustar,  iw,  ix,  upfactor,  iflag, 

&  sigmaudp,  sigmavdp,  sigmawdp,  Tlw,  sumdw,  dnfactor,  dwfactor) 
c  this  subroutine  computes  particle  velocity  and  position  for  different 
c  values  of  Ri#,  etc. 
c 

REAL  Km,  L,  imdamu,  lmdamv,  lmdamw 
DOUBLE  PRECISION  Linv,  ruu,  rvv,  rww,  pu,  pv,  pw 
c  Get  mesoscale  and  boundary  layer  flow  parameters 
c  Initialize  variables 

zi  =  2000  -  6000* (Linv  +  .2) 

zovrzi  =  z/zi 

zovrL  =  z*Linv 

ziovrL  =  zi*Linv 

abszovrL  =  abs( zovrL) 

L  =  l./Linv 

CALL  MESOFLOW (  z,  zi,  zO,  Linv,  windspd2,  ustar,  Km,  Ri, 

&  u,  v,  w,  Theta2,  dudz,  dvdz,  wk,  D,  R,  wptpZ,  BVF,  iw, 

&  BLS,  sigmaU2,  sigmaV2,  sigmaW2,  tke,  Wz,  i,  wk,  ustar. 
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&  bdthetdz, Thstar,  Deltheta) 

d  IF ( i . It . iw) WRI^E ( 6 , * ) '  zov-zi,  zovrL,  ziovrL  '  ,  zovrzi ,  zovrL,  ziovrL 

c  Get  A,  lambdas,  ou,v,wdp,  and  ouoyancy  length  scale,  BLS 
IF  (  zovrL  .It.  0.0  )  THEN 

sigmaudp  =  ustar*(12.0  +  0. 5*zi/abrM'L)  )  **0. 3333333 
IF(i. It. iw) WRITE  ( 6 , * ) ’ ou  "  =u* ( 12  + . azi/ | L | ) **1/3  , 

&  sigmaudp, ustar , zi, L 

IF  (  abszovrL  .le.  1.  ) 

&  A  =  0 . 31* ( 1 .  -  3 . * zovrL) **{-0.333333)*(l.  -  15 . *zovrL) **0 . 25* 
Si  (0.55  +  0.38*zovrL  ) 

IF  (  0. l*abs( ziovrL)  .gt.  abszovrL  .and.  abszovrL  .gt.  1.) 

Si  A  =  0 . 05*  ( 1 .  -  3 .  *  zovrL)  **(-0.333333)*(l.  -  15  .  *zovrL)  **0 . 25 
A  =  MAX  (  A,  0.07) 

A  =  MIN  (  A,  0.18) 
lmdamu  =  1.5*zi 

IF  (  z  .le.  0.1*zi  .and.  abs(zovrl)  .gt.  1.)  Imdamw  =  5.9*z 
IF  (  z  .le.  abs(L)  )  Imdamw  =  j/(  0.55  +  0.38*zovrL  ) 

IF  (  0.1*zi  .It.  z  .and.  z  .It.  zi  ) 

&  Imdamw  =  1.8*zi*(  1.  -  exp(-4. 0*zovi ri) 

&  -  0.0003*exp(  8.0*zovrzi  )  ) 

sigmawdp  =  Km/ (A*lmdamw) 

IF(i.lt.iw) 

&  WRITE ( 6, * ) ’ at  z , ow *' =Km/ (A*lmdamw) z, sigmawdp, Km, A, Imdamw 

ELS'  *F  (  zovrL  .ge.  0.)  THEN 
sigmaudp  =  2.3*ustar 
lmdamu  =  0 . 7*sqrt ( zovrzi) *zi 
es  =  70.0 

IF  (z. It. 205.0)  es  =  0. 35*z 
Imdamw  =  MAX (  z,  2.9*es) 

Ric  =  0.25 

Rifactor  =  ((Ric  -  1 . 25*Ri ) /Ric ) **0 . 58 

IF(i. It. iw) WRITE (6, *) 1 Rif= ( (Ric-Ri ) /Ric ) ** . 58 ’ , Rifactor, Ric, Ri 
shear  =  sqrt(dudz**2  +  dvdz**2) 
sigmawdp  =  1 . 2*es*Rifactor*shear 
IF ( i . It . iw ) WRITE ( 6 , * ) ’ow' ' =1 . 2esRi , shear ’ , 

Si  sigmawdp,  es,  Rifactor,  shear 

BLS  =  sigmawdp/ BVF 
ENDIF 

sigmavdp  =  sigmaudp 
IF  ( ikamada  .eq.  1)  THEN 

Use  values  derived  from  Sorbjan  and  Kamada  from  mesoscale  subroutine 
IF(i.lt.iw)WRITE(6,*)'ou2,av2,aw2’, sigmaU2 , sigmaV2 , sigmaV72 
sigmaudp  =  sqrt ( sigmaU2 ) 
sigmavdp  =  sqrt ( sigmaV2 ) 
sigmawdp  =  sqrt ( sigmaW2 ) 

IF ( i . It . iw) WRITE ( 6 , * ) ' iKmda  ou,v,w’ * ' , sigmaudp, sigmavdp, s : gmawdp 
To  get  vertical  integral  length,  see  pps.  244-247  Lectures  on  Air 
Pollution 

aln  =  0.41*z 

Imdamw  =  0 . 27*sigmawdp/BVF 
Imdamw  =  l./(l./aln  +  1. /Imdamw) 

IF ( i . It . iw) WRITE ( 6 , * ) ' iKmda  aln, Imdamw' , aln, lmda . w 
ENDIF 

IF  (noskew  .eq.  1)  tke  =  0.5*(sigmaU2  +  sigmaV2  +  sigmaW2) 
lmdamv  =  lmdamu 

IF(i.lt.iw)WRITE(6,*) ' lmdamu , lmdamv, Imdamw ' , lmdamu , lmdamv, Imdamw 
Get  mean  velocity 

V  =  sqrt(  u*u  r  v*v  +  w*w  ) 

Get  betas 

IF(i.lt.iw)  WRITE  ( 6 ,*)’ u , v , w, sigmaudp ', u , v,w, sigmaudp 
IF(i.lt.iw)  WRITE  ( 6 ,*)' s igmavdp , sigmawdp ', s igmavdp , s igmawdp 
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betau  =  0. 6*V/sigmaudp 
betav  =  0. 6*V/sigmavdp 
betaw  =  0. 6*V/sigmawdp 

IF{ i. It . iw) WRITE ( 6, * ) ’at  z, flw= . 6V/cw* ' ' , z, betaw, V, sigmawdp 
Get  Lagrangian  time  scales 

IF( i. It . iw) WRITE ( 6, * ) ' betau, lmdamu, V' , betau, lmdamu, V 
Tlu  =  0.2*betau*lmdamu/V 
Tlv  =  0.2*betav*lmdamv/V 
Tlw  =  0.2*betaw*lmdamw/V 

IF( i. It . iw) WRITE ( 6, * ) ' Tlw= . 2*flw*lmdamw/V' , z, Tlw, bet aw, lmdamw, V 
IF  (  ikamada  .eg.  1)  THEN 

IF  (  Linv  .ge.  0.  )  Tlw  =  lmdamw/ sigmawdp 
IF  (  Linv  .It.  0.  )  THEN 
Tlw  =  0.3*zi/wk 
AspectR  =  2.0  -  40.0*Linv 
Tlu  =  AspectR*Tlw 
Tlv  =  Tlu 
END  IF 

IF( i. It . iw) WRITE ( 6, * ) ' Tlw= lmdamw/ cw ' ' ’ ,Tlw, lmdamw, sigmawdp 
ENDIF 

IF  (  Tlu  .It.  0.02*dt  )  Tlu  =  0. 02*dt 
IF  (  Tlv  .It.  0.02*dt  )  Tlv  =  0.02*dt 
IF  (  Tlw  .It.  0.02*dt  )  Tlw  =  0.02*dt 
Get  Autocorrelation  function 

IF ( i . It . iw) WRITE ( 6, * ) 1 dt , TLu, TLv, TLw 1 ,  dt,  Tlu,Tlv,TLw 
Ru  =  exp(  -dt/Tlu) 

Rv  =  exp(  -dt/Tlv) 

Rw  =  exp(  -dt/Tlw) 

IF(i. It. iw) WRITE  ( 6 , * ) ’at  z,Rw=exp( -dt/Tlw) • ,z,Rw,dt,Tlw 
IF( i. It . iw) WRITE  ( 6, * ) ’ Ru,Rv,Rw, ikamada ’ , Ru,Rv,Rw, ikamada 
Get  standard  deviations  of  velocity 

sigmautp  =  sigmaudp*sqrt (  (1.  -  Ru**2)  ) 
sigmavtp  =  sigmavdp*sqrt (  (1.  -  Rv**2)  ) 
sigmawtp  =  sigmawdp* sqrt (  (1.  -  Rw**2)  ) 

IF(i.lt.iw) 

&  WRITE ( 6, * ) ’ at  z, ow' ' '=ow' '/ ( l-RwJ )’, z, sigmawtp, sigmawdp, Rw 

IF(i.lt.iw) 

&WRITE (6,*),CTu,'’,av',',ow',',, sigmautp, sigmavtp, sigmawtp 
Get  random  normal  fluctuation  velocities,  utp  &  vtp,  at  time,  t-dt 
ruu  =  abs(pu/3.d0) 
rw  =  abs(pv/3.d0) 
rww  =  abs(pw/3.d0) 

IF(i.lt.iw)  WRITE  (6, *) ' ruu , rw , rww ’ ,ruu,rvv,rww 
call  NORNG(  ruu,  pu  ) 
call  NORNG(  rw,  pv  ) 

IF (  pw  .ge.  0.0  ) wtp  =  upfactor* (pw* ( sigmaw2  -sigmawtp)  +  wbuoy) 
utp  =  sigmautp*pu 
vtp  =  sigmavtp*pv 

IF  (  noskew  .ne.  1  .and.  Linv  .It.  0.  )  GOTO  250 
call  NORNG (  rww,  pw  ) 
wtp  =  sigmawtp*pw  +  wbuoy 
IF(i.lt.iw)  WRITE  ( 6 , * ) ' pu, pv, pw ' , pu , pv, pw 
IF(i.lt.iw)  WRITE  ( 6 , * ) ' utp, vtp, wtp ' , utp, vtp, wtp 
IF(i.lt.iw)  WRITE  (6,*)'w'*'  =  cw ''' *pw ’, wtp, sigmawtp, pw 
GOTO  500 
50  CONTINUE 

IF(i.lt.iw)  WRITE  ( 6 , * ) ’ pu , pv, pw ’ , pu , pv, pw 
IF(i.lt.iw)  WRITE  ( 6 , * ) ’ utp, vtp, wtp ’ , utp, vtp, wtp 
if  (  iMcNid  .ne.  1  )  GOTO  400 
c  modify  wtp  for  skewness  according  to  McNider  method 
iwp  =  0 
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iwm  =  0 

300  CONTINUE 

c  Get  random  normal  fluctuation  velocity,  wtp,  at  time,  t-dt 
call  NORNG (  rww,  pw  ) 
wtp  =  sigmawtp*pw 


IF 

( 

wtp  .ge.  0. 

.  and. 

iwp  .eq. 

0 

) 

wp 

=  wtp 

IF 

{ 

wtp  .ge.  0. 

)  iwp 

=  1 

IF 

( 

wtp  .le.  0. 

.  and. 

iwm  .eq. 

0 

) 

wm 

=  wtp 

IF 

( 

wtp  .le.  0. 

)  iwm 

=  1 

IF 

( 

iwp  .eq.  0 

.or.  iwm  .eq.  0 

) 

GOTO 

300 

c  Get  S  function 

IF  (  zovrL  .gt.  0.  )S  =  0.1  -  0.2*zovrL**0.2 

IF  (  zovrL  .le.  0.  )S  =  0.1  +  (0.6/(  0. 68* ( 1 . -15 . *zovrL) ** (-0. 25 ) 

&  -  1.8*zovrL  )  ) 

d  IF(i.lt.iw)  WRITE  ( 6, * ) ' iwp, iwm, S ' , iwp, iwm, S 

c  Get  alpha  and  eta 

alpha  =  -  0.028  -  0.6*abs(S) 
eta  =  0.54*abs(S) 

c  Get  wtp(t-dt)  111  Pielke  prints  this  on  p.  178  as  wdp(t-dt)  CHECK  THIS 
c  by  looking  at  McNider  paper!!!!! 

IF  (  zovrL  .le.  0.  )  wtp  =  alpha*wp  +  wm/alpha  -  eta 

IF  (  zovrl  .gt.  0.  )  wtp  =  wp/alpha  +  alpha*wm  +  eta 

d  IF(i.lt.iw)  WRITE  ( 6 ,*)' alpha, eta, wtp ’, alpha, eta, wtp 

GOTO  500 
400  CONTINUE 

IF  (  iKamada  .ne.  1  .and.  Linv  .It.  0.  )  GOTO  500 
c 

c  Try  Kamada-Berkowicz  double  gaussian  skewness  approx  instead  of  McNider. 
c  This  approach  assumes  that  mean  absolute  vertical  fluctuation  velocity 
c  is  0 . 55wk,  vertically  averaged  over  the  BL  depth,  where  wk  generalizes 
c  w*  to  mixed  forced/free  convection.  We  let  the  up/downdraft  volume 
c  ratio  in  the  convective  BL  "  0.4/0. 6,  with  mean  updraft  velocity 
c  *  1.08ow  and  the  mean  downdraft  velocity  "  0.87ow.  We  assume  that  the 
c  two  gaussian  velocity  distributions  are  displaced  from  zero  by  this 
c  amount  with  60%  of  the  distribution  in  the  downdraft  volume  and  40% 
c  in  the  updraft  volume.  We  assume  that  departures  from  Gaussian  are 
c  small  for  stable  cases.  (See  pps.  199-201,  Ch.  4  Lectures  on  Air 
c  Pollution  Modeling) 
c 

utp  =  pu*sigmaudp*sqrt ( 1  -  Ru**2) 
vtp  =  pv*sigmavdp*sqrt ( 1  -  Rv**2) 
wbuoy  =  Rw*wbuoy  -  ( 9 . 801/Theta2 ) *Deltheta*dt 
sumwbuoy  =  sumwbuoy  +  wbuoy 
call  NORNG (  rww,  pw  ) 

IF  (  pw  .ge.  0.0  )  ikount  =  ikount  +  1 
d  IF(i.lt.iw)  WRITE  ( 6 ,*)' ikount , icount ', ikount , icount 

IF  (  ikount  .eq.  icount  )  THEN 
pw  =  -pw 
ikount  =  0 


END  IF 

The  above  sets 


c  E. 

g. ,  if 

icount 

420 

CONTINUE 

c 

IF  ( 

pw  . ge . 

c 

IF  ( 

pw  .le. 

IF  { 

pw  . ge . 

& 

wtp  = 

IF( 

pw  .le. 

& 

wtp  = 

c 

IF  ( 

pw  .le. 

an  (icount  +  1) /(2*icount)  chance  of  updraft  &  an 
(icount  -  1 ) / ( 2*icount )  chance  of  a  downdraft. 
5,  then  updraft /downdraft  chances  are  60%/40%. 


0.0  )wtp  =  upf actor*sigmawtp*pw  +  wbuoy 
0.0  )wtp  =  dnf actor*sigmawtp*pw  +  wbuoy 
0.0  ) 

wtp  =  upf actor*pw*sigmawdp*sqrt ( 1  -  Rw**2)  +  wbuoy 

0.0  ) 

wtp  =  dnfactor*pw*sigmawdp*sqrt ( 1  -  Rw**2)  +  wbuoy 
0.0  )  wtp  =  upf actor*pw*sigmawdp  +  wbuoy 
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c  IF(  pw  .le.  0.0  )  wtp  =  dnf actor* pw*sigmawdp  +  wbuoy 

sumw  =  sumw  +  wtp 

d  IF( i. It . iw) WRITE ( 6, * ) ' wtp, pw' , wtp, pw 

d  IF( iflag  .eq.  iw) 

d  &  WRITE( 6, * ) ' wbuoy, Rw, 6©, sumw' ,  wbuoy, Rw,Deltheta, sumw 

d  IF(iflag  .eq.  ix)  WRITE ( 6 ,*)’ sumwbuoy, sumw' ,  sumwbuoy, sumw 

500  CONTINUE 

c 

c  Get  fluctuation  velocities 
udp  =  udpold*Ru  +  utp 
vdp  =  vdpold*Rv  +  vtp 
wdp  =  wdpold*Rw  +  wtp 
IF  (ikamada  .eq.  1)  THEN 

udp  =  udpold  +  utp  -  (1.  -  Ru)*udpold 

vdp  =  vdpold  +  vtp  -  (1.  -  Rv)*vdpold 

wdp  =  wdpold  +  wtp  -  (1.  -  Rw)*wdpold 

ENDIF 

d  IF(i.lt.iw) 

d  &  WRITE  (6,*) 'at  z, wdp=wdpold*Rw+wtp ’ , z, wdp, wdpold, Rw, wtp 
c  Get  position 

x  =  x  +  udp*dt 
y  =  y  +  vdp*dt 
z  =  z  +  wdp*dt 
IF  (  z  .It.  10*z0  )  THEN 
z  =  abs(z)  +  10.*z0 
wdp  =  abs(wdp) 

ELSEIF  (  z  .gt.  zi  )  THEN 
z  =  2.0*zi  -  z 
z  =  MAX ( z ,  0.9*zi) 
wdp  =  -abs(wdp) 

ENDIF 

udpold  =  udp 

vdpold  *  vdp 

wdpold  =  wdp 

iflag  =  iflag  +  1 

IF  (iflag  .gt.  ix)  THEN 

d  WRITE  (6, 1110)time,Iiinv,Ri,udp,  vdp,wdp,x,y,  z 

d  WRITE  (6,1120)Km,  V,  ustar,  wpTpZ,  sigmaW2,  tke,  Tlw,  Rw,  Wz 

iflag  =  1 
ENDIF 

*********************************************************************** 

c  debug  formats 

dlOOO  FORMAT  (lOx,  2gll.3) 

dlOlO  FORMAT  (lOx,  3fll.3) 

dl020  FORMAT  (lOx,  4fll.3)  ■ 

dl030  FORMAT  (lx,  i5  ) 

dl040  FORMAT  (lOx,  3gll.3,  i5) 

dl050  FORMAT  (lOx,  4fll.3,  i5) 

dl055  FORMAT  (lOx,  5fll.3,  i5) 

dl060  CONTINUE 

dlllO  FORMAT  (lx, ' tm  1/L  Ri  u  v  w  x  y  z  ’ , f 6 . 1, 2f 6. 2 , 3f 6 . 1, 3f 7 . 1 ) 
dll20  FORMAT  (lx, 'Km  V  u*  Q  owJ  e  Tw  Rw  Wz  ' , f 6 . 1 , 5f 6 . 2 , f 6 . 1 , 2f 6 . 2 ) 
RETURN 
END 

it********************************************************************* 

SUBROUTINE  MESOFLOW  (  z,  zi,  zO,  Linv,  windspd2,  ustar.  Km,  Ri, 

&  u,  v,  w,  Theta2,  dUtotdz,  dvdz,  wk,  D,  R,  wpTpZ,  BVF,  iw, 

&  BLS,  sigmaU2,  sigmaV2,  sigmaW2,  tke,  Wz,i,wk,  ustar, 

&  bdthetdz,  Thstar,  Deltheta) 

c  This  subroutine  computes  the  mesoscale  flow  parameters  used  as 
c  inputs  to  the  standard  mcnider  particle  calculations.  It  also  computes 
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some  alternate  values  for  the  ou,  v,  w,  and  other  quantities. 

REAL  k,  LI,  InzovzO,  LovZi,  Km,  L,  lepsilon,  lk 
DOUBLE  PRECISION  Linv 

DATA  g,  k,  zl,  f  /  9.801,  0.4,  2.0,  r.0001  / 

IF ( i . It . iw) WRITE { 6 , * ) 'in  MESOFLOW  i  =  ',i 
Zero  out  some  variables 
u  =  0 . 
v  =  0. 
w  =  0. 
dUdZ  =  0. 
dVdZ  *  0. 
wstar  =  0.0 
Wk  =  0. 

Ri  =  0. 

Km  =  0. 
tke  =  0. 

BVF  =  0. 

•  BLS  «  0. 

Initialize  some  values 
D  =  0.1 
R  =  0.2 
zOinv  =  l./zO 
InzovzO  =  alog( z*z0inv) 

IF  (  z  .It.  zO  )  z  =  z  +  zO 

zinv  =  l./z 

ziinv  =  l./zi 

L  =  l./Linv 

gzi  =  9.801*zi 

zlovL  =  zl*Linv 

Theta2  =  290.0 

zovL  =  z*Linv 

Get  psim  &  psih  for  surface  layer.  If  stable, 
psim  =  -5 . *zlovL 
psih  =  psim 

IF( i. It . iw) WRITE ( 6, * ) ' psim=psih=-5zl/L • , psim, psih, zlovL 
If  unstable, 

IF  (  Linv  .It.  0.)  THEN 

phim  =  (1.  -  28. *zlovL) ** (-0. 25 ) 

IF(  i.  It .  iw)  WRITE ( 6,  * )  '  phim=( l-28zl/L)  **  (-**)  '  ,  phim,  zlovL 
Get  psim  by  Kamada  (unpub) 

IF  (  zlovL  .It.  0  .and.  zlovL  .gt.  -0.01)  zlovL  =  -0.01 
IF( i. It. iw)WRITE(6, *) ' just  before  psim  =  ' 

psim  =  (.0159651383  -  5 . 4151107*zlovL) / 

&  (1.  -  3 . 590022 31* zlovL  -0 . 799168457*zlovL* *2 ) 

IF ( i. It . iw) WRITE ( 6, * ) ' psim= ( a-czl/L) / ( l-bzl/L-d( zl/L) 2 ' , psim, zlovL 
Also  from  Dyer  and  Bradley  (1982)  BLM  22,  3-19,  we  have  for  psi 

IF( i. It . iw) WRITE ( 6, * ) ' just  before  psih  =  ’ 

psih  =  2.*alog(  0.5*(1.  +  sqrtf  1.  -  14.*zlovL  ))  ) 

IF(  i .  It .  iw)  WRITE  ( 6,  *  )  '  psih=2Ln(lj  ( 1+/  ( 1-14 zl/L)  )  )  '  ,psih,  zlovL 
ENDIF 

Get  windspeed  drag  coefficient  and  friction  velocity 
IF( i. It . iw) WRITE ( 6, * ) ' just  before  sqrtCdn  =  ' 
sqrtCdn  =  0.4/(  alog( zl*z0inv)  -  psim) 

IF( i. It . iw) WRITE ( 6 , * ) '/ Cdn=. 4/ (Ln ( zl/zO ) -psim) ' , sqrtCdn, zl , zO, psim 
ustar  =  MAX (  sqrtCdn*windspd2 ,  0.01  ) 

IF( i . It . iw) WRITE ( 6, * ) ' u*=MAX(/ Cdn*V2 ,0.01) ' , ustar, sqrtCdn, windspd2 

For  dispersion,  get  vertical  profiles  of  relevant  variables,  using 
similarity  theory  from  (Sorbjan  , Structure  of  the  Atmos  BL,  Ch.  4)  and 
(Panofsky  and  Dutton,  Atmos  Turbulence).  First  get  oft  used  variables 


145 


zovzi  =  z*ziinv 
zOovL  =  zO*Linv 
ziovL  =  zi*Linv 
Lovzi  =  L*ziinv 
ustar2  =  ustar*ustar 
ustar3  =  ustar*ustar2 
wpTpO  =  -0.25508*Linv*ustar3*Theta2 
d  IF( i. It . iw) WRITE  (6,*)'w0=  =-. 255u**30/L wpTpO , ustar, Theta2 , L 

Tatar  =  -wpTpO/ustar 
IF  (  Linv  .It.  0.)  THEN 

Wstar  =  (gzi*wpTpO/Theta2 )**. 33333 

d  IF ( i . It . iw) WRITE ( 6 , * ) ' w*= (gziwe/0) **1/3 ' , wstar , gzi, wpTpO, Theta2 

c  Replace  w*  with  wk  for  mixed  forced/free  convection  (Kamada, 

c  NPS-PH-92-007 .  Here,  add  only  surface  layer  shear  production. 

Ar  =  (1.  -  150. 0*z0ovL)**0. 333333  -  1 
d  IF( i. It. iw) WRITE  (6,*)'Ar=  (1  -  150z0/L) ** ( 1/3 )  -  1* 

d  IF( i. It. iw) WRITE  (6,*)Ar,  zOovL 

d  IF( i. It . iw)WRITE( 6, * ) ’ just  before  Wk  =  ' 

Wk  *  wstar*(l.  -  (  3.125  -  2.5*alog(Ar)  ) *Lovzi  )**0. 33333 
d  IF(i. It. iw) WRITE  < 6 , * ) ’ Wk=w* ( 1-3 . 124-2 . 51n ( Ar ) )L/zi) ** ( 1/3 )  ’ 

d  IF ( i . It . iw) WRITE  ( 6 ,* )wk, wstar , Ar, Lovzi 

Thstar  =  -wpTpO/wk 

c  Decide  what  algos  to  use,  based  on  z/L. 

END  IF 

IF  (  zovL  .It.  -0.5  )  icase  =  1 

IF  (  zovL  .It.  0.0  .and.  zovL  .ge.  -0.5  )  icase  =  2 

IF  (  zovL  .ge.  0.01  )  icase  =  3 

GOTO  (  330,  400,  430  )  icase 
330  CONTINUE 

********* 

c  Get  Km,  &  ou,v,w  for  forced/free  convective  BL.  Ri  and  shear  not 
needed. 

********* 

c  Get  potential  temperature  at  height  z  above  -L/2 

R  =  abs(R) 

zovzil3  =  zovzi** . 33333 
zovzi23  =  zovzil3*zovzil3 
tpl  =  1.0  -  zovzi 
tp4  =  tpl  +  D 
tp413  =  tp4**. 33333 
tp423  =  tp4**. 66667 
Thstar  =  -wpTpO/wk 
tp23  =  tpl**. 666667 
R23  =  R**. 666667 
Oldtheta  =  Theta 

Theta  =  Theta2  +  Thstar* (  tp23/zovzil3  +  R23*zovzi23/tp413  ) 

IF  (  i.eq.  1)  Oldtheta  =  Theta 
Deltheta  =  Theta  -  Oldtheta 

d  IF( i. It. iw)WRITE  ( 6 , * ) ' SO, O, Old© ' ,  Deltheta,  Theta,  Oldtheta 

c  temperature  and  buoyancy  fluxes  at  height,  z 

wpTpZ  =  wpTpO* (1.  -  1.2* zovzi) 

B  =  g/Theta2 
BwpTpZ  =  B*wpTpZ 

c  Get  vertical  velocity  variance  at  height  z  above  -L/2 

sigmaW2  =  wk**2* ( 1. l*zovzi23*tp23  +  R23*zovzi23*tp423 ) 
d  IF( i. It. iw) WRITE  ( 6 , * ) ' ow» =wk J ( 1 . 1 ( z/zi ) **2 /3 ( 1-z/ zi ) **2 /3  +  ’ 

d  IF< i. It. iw) WRITE  (6,*)’R**(2/3)(z/zi)**(2/3) ( 1-z/zi+D) ** ( 1/3  )  ' 

d  IF ( i . It . iw) WRITE  ( 6 , * ) sigmaW2 , wk , zovzi , R, D 

sigmaW  =  sqrt ( sigmaW2 ) 

c  Get  TKE  dissipation  rate  at  height  z  above  -L/2,  using  similarity 

tmplO  =  Wk**3*ziinv 
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epsilonO=  0.75*tmpl0 
phiEp  =  0.75*tpl  +  R*zovzi 
epsilon  =  phiEp*tmplO 

IF (i. It. iw) WRITE  (6,*)'e  =  *e*wk**3/zi’ , epsilon, phiEp,wk, zi 
Get  au1 ,av2 

wpsqovE  =  0.333*(  (2.*epsilon0  -  epsilon) /epsilonO  + 

&  BwpTpZ /epsilon  ) 

IF ( i . It . iw ) WRITE ( 6 , * ) *w' »/e=(2eO-e)/3eO  +Gw'0,z/e' , wpsqovE 
IF( i. It . iw)WRITE ( 6, * ) ' cO,  fiw'0'z*,  epsilonO,  BwpTpZ 
sigmaU2  =  0. 5555*sigmaW2* (  2. /wpsqovE  -  1.) 

c  sigmaU2  =  MAX(  sigmaU2,  .33333  ) 

sigmaV2  =  0.8*sigmaU2 

d  IF( i. It. iw)WRITE  ( 6 , * ) sigmaU2 , wk, zovzi , sigmW2 

c  The  following  tke  parameterization  is  OK  for  low  shear 

tke  =  0.5*(  sigmaU2  +  sigmaV2  +  sigmaW2  ) 
d  IF( i. It . iw) WRITE  (6,*)’e  =  ow2/2  +•  3oul/2 ’ ,tke, sigmaW2, sigmaU2 

tkel2inv  =  l./sqrt(tke) 
tke32inv  =  tkel2inv**3 

c  Get  dissipation  &  mixing  length  scales  &  vertical  diffusion 
c  coefficient 

lepsilon  =  l./(  zinv  +  1 . 4*epsilon*tke32inv  ) 
d  IF{ i. It . iw) WRITE ( 6, * ) • le  =  l/(l/z  +1.4e/e  ) 1 ,  lepsilon, z, epsilon, tke 
lk  =  MIN  (  z,  3 . 0*lepsilon*sigmaW2/tke  ) 
d  IF( i. It . iw) WRITE  (6,*)'lk  =  MIN(3l6*ow2/tke,  z) ', lk, lepsilon, tke, z 

c  The  mixing  length  scale  and  diffusion  coefficients  are  Kamada 

c  variants 

Km  =  0.44*lk/tkel2inv 

d  IF( i. It. iw)WRITE  (6,*) ’ Km  =  0 . 44*lk*sqrt (tke) ’ ,Km,lk,tke 

c  Get  temperature  and  buoyancy  gradients 

c  dthetadz  =  . 6*Thstar*ziinv* (tp23/zovzi23**2  -2 . *R23*zovzi23/tp423 ) 

dthetadz  =  .6*Thstar*ziinv*(tp23/zovzil3  -2 . *R23*zovzi23/tp423 ) 
Bdthetdz  =  B*dthetadz 
GOTO  460 
400  CONTINUE 

c  get  Km  for  unstable  surface  layer.  Shear  and  Ri  not  needed  for 
unstable. 

Km  =  k*ustar*z/phim 
tkeusl  =  ( 5 . 6*Km*zinv) **2 
Ri  =  zovL 

phil  =  (  12.  -  0. 5*ziovL  >**.333333 
phi2  =  phil 

phi3  =  {  1.  -  14.*zovL  )**.25 
phiEp=  (1.  +  .75* ( -zovL) **.6667)**1.5 
epsilonO  =  ustar3*phiEp*zinv 
GOTO  440 
430  CONTINUE 
********** 

c  get  Km,  Ri,  and  other  parameters  for  neutral  to  stable  BL 
********** 

h  =  z  i 

IF  (  zi  .It.  -98.  )  h  =  sqrt (k*ustar*Linv/f ) 
d  IF( i. It . iw) WRITE ( 6, * ) ' h=zi,or  / ( ku*/ ( fL) ) ' , h, zi, f 

hinv  =  l./h 
zovH  =  z*hinv 
zsqovHL  =  zovL*zovH 
zsqovH2  =  zovH*zovH 
alpha  =  2.0  -  10.0*Linv 
beta  =  3.0  -  20.0*Linv 

d  IF(  i.  It.  iw)WRITE(6,  *  )  'a, 13,  z/H'  ,  alpha,  beta,  zovH 
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tp8  =  1.  -  zovH 
tp8alpha  =  tp8**alpha 

d  IF(  i.  It .  iw)WRITE  ( 6,  *  )  '  ( 1-z/H) ,  ( 1-z/H)  **a  '  ,  tp8ftp8alpha 

tp8al2  =  sqrt (tp8alpha) 
tp8beta  =  tp8**beta 

d  IF ( i . It . iw) WRITE ( 6 ,  *  ) * ( 1-z/H) **a/2 , (l-z/H)**Bl , tp8al2 , tp8beta 

phil  =  2.2  *tp8al2 
phi2  =  phil 
phi3  =  1.6*tp8al2 
c  local  friction  velocity 

ul  =  ustar*sqrt (tp8alpha) 

d  IF( i. It . iw ) WRITE ( 6 , * ) ' ul=u* (l-z/H)**(a/2) ' , ul, ustar , zovH 

c  local  Obukhov  length 

LI  =  L*tp8** ( 1 . 5*alpha  -  bef\) 

IF(  i .  It .  iw)  WRITE  ( 6 ,  *  )  '  L1=L  ( 1-z/H)  **(  3a/2-!i)  '  ,  LI  ,L,  zovH,  alpha, 
&  beta 


total  shear 

dUtotdz  =  4.7*ul/(k*Ll) 

dUtotdz  =  2 . 5*ustar*zinv* ( 1 . +4 . 7*zovL) ** (0. 5*alpha) *tp8al2 
IF  (  i .  It .  iw) WRITE  ( 6 ,  * )  '  dU/dz=4 . 7ul/ ( kLl )  ’  , dUtotdz , ul , k, LI 
temperature  and  buoyancy  fluxes  at  height,  z 
wpTpZ  =  wpTpO*(l.  -  zovzi)**beta 
Brunt  Vaisala  frequency  at  z 
tp9  =  (1.  +  3 . 7*zovL) *zinv 

IF( i . It . iw) WRITE ( 6, * ) *  tp9={ 1+3 . 7z/L) /z ‘ , tp9, zovL, zinv 
BVF  =  4 . 3*ustar*Tstar**2*tp9*tp8beta*tp8alpha 

IF(i. It. iw)WRITE(6, *) ' BVF=4 . 3u*Q*  2 ( 1+3 . 7z/L) (1-z/H) **(a+B) /z1 
&  BVF, ustar, Tstar, zovH, zovL, z 
TKE  dissipation  rate  at  z 

phiEp  =  3 . 6*tp9*tp8beta*zinv 
epsilon  =  phiEp*ustar3 
Richardson  number  at  z 

tplO  =  l./(l.  +  5.*zovL  ) 

Ri  =  zovL*tplO 

momentum/heat  dif fusivities  at  z 
Km  =  k*ustar*z*tp8*tpl0 

IF ( i . It . iw) WRITE ( 6, * ) ' Km=ku*z ( 1-z/H) /( 1+5 z/L) ' , Km, k, ustar , z. 


&  zovH, zovL 

**************************** 


0  CONTINUE 

Get  ou,v,w,  tke,  total  horizontal  velocity,  &  temperature 
sigmaU  =  ustar*phil 
sigmaU2=  sigmaU*sigmaU 
sigmaV  =  0.894*sigmaU 
sigmaV2=  0.5*sigmaU2 
sigmaW  =  ustar*phi3 
sigmaW2  =  sigmaW**2 

tke  =  0.5*(sigmaU2  +  sigmaV2  +  sigmaW2) 

IF ( i. It . iw) WRITE ( 6, * ) ' tke=oui2/2 ' , tke, sigmaU2 , sigmaV2 , sigmaW2 
IF  (  icase  .eq.  2  )  tke  =  tkeusl 
Oldtheta  =  Theta 

Theta  =  Theta2  +  2.5*Tstar*(  InzovzO  -  psih  ) 

IF  (  i.  eq.  1)  Oldtheta  =  Theta 
Deltheta  =  Theta  -  Oldtheta 

IF(  i.  It.  iw)WRITE  ( 6 ,  * )  '  6<j>,  0,  Old© '  ,  Deltheta,  Theta,  Oldtheta 
i0  CONTINUE 

Vtotal  =  2 . 5*ustar* ( InzovzO  -  psim) 

IF ( i . It . iw) WRITE ( 6 , * ) ' V=2 . 5u* (Ln ( z/zO ) -psim) ) ' , Vtotal , ustar , 

&  lnzovzO,psim 

u  =  Vtotal 

d  IF ( i . It . iw) WRITE ( 6, * ) ' 0=02+2 . 50* ( Ln ( z/zO ) -psih) ' , theta , theta2 , 
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d  &  tstar, lnzovzO,psih 

********************************************************************** 

c  debug  formats 
dlOOO  FORMAT  (lOx,  2gll.3) 
dlOlO  FORMAT  ( lOx,  3fll.3) 
dl020  FORMAT  (lOx,  4fll.3) 
dl030  FORMAT  (lx,  i5  ) 
dl040  FORMAT  (lOx,  3gll.3,  i5) 
dlOSO  FORMAT  (lOx,  4fll.3,  i5) 
dl055  FORMAT  (lOx,  5fll.3,  iS) 
dl060  CONTINUE 
RETURN 
END 

********************************************************************** 
SUBROUTINE  NORNG ( r , p ) 
c 

c  This  subroutine  generates  a  sequence  of  numbers  normally  and  randomly 
c  distributed  over  the  interval  -3  to  3  from  uniformly  distributed  random 
c  numbers,  by  the  method  of  linear  approximation  to  the  inverse  of  the 
c  accumulative  normal  distribution  function 
c 

DOUBLE  PRECISION  r,p,  y<6),  x(6),  3(5) 

DATA  y/O.dO,  0.0228d0,  0.0668d0,  0.1357d0,  0.2743d0,  0.5d0  /, 

&  x/  -3 . OldO ,  -2 . OdO,  -1.5d0,  -l.OdO,  -0.6d0,  O.OdO  /, 

&  s/  43 . 8596d0 ,  11.3636d0,  7.25689d0,  2.891352d0,  2.65887d0  / 

CALL  STRNUM(r) 
p  =  r 
i  =  1 

IF  (  p  .gt.  0 . 5d0  )  p  =  l.OdO  -  r2  IF  (  p  .It.  y(i+l)  )  GOTO  8 

i  =  i  +  1 
GOTO  2 

8  p  =  (  (  p  -  y(i)  )*3(i)  +  x(i)  ) 

IF  (  r  .ge.  0.5d0  )  p  =  -p 
cd  write  (6,*) 'random  numbers,  r  &  p',r,p 
RETURN 
END 

************************************************************************ 
SUBROUTINE  STRNUM(r) 

c 

c  This  subroutine  generates  a  sequence  of  numbers  which  are  randomly  and 
c  uniformly  distributed  over  the  unit  interval, 
c 

DOUBLE  PRECISION  pi,  r 
bb  =  l.dO 
pi  =  r*317.d0 
c  WRITE  ( 6 , * ) ' pi ' , pi 

r  =  abs(  pi  -  (  INT(pl/bb) *bb)  ) 

RETURN 

END 

************************************************************************ 
FUNCTION  MAX  (a,  b) 
max  =  b 

IF  (  a  .gt.  b)  max  =  a 
END 

************************************************************************ 
FUNCTION  MIN  (a,  b) 
min  =  b 

IF  (  a  . It .  b)  max  =  a 
END 

************************************************************************ 
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